1、Eigen每个头文件的作用
2、Eigen使用基础
2.1 Hello World
#include<iostream> #include<Eigen/Dense> // Eigen头文件,包含Eigen库里面所有的函数和类 int main() { Eigen::MatrixXd m(2,3); // MatrixXd:动态数组,初始化时指定数组的行列 m(0,0) = 3; m(1,0) = 2.5; ...; std::cout << m << std::endl; // Eigen重载了<<运算符,可以直接输出Eigen矩阵的值 }
2.2 Matrices and Vectors
#include<iostream> #include<Eigen/Dense> using namespace std; using namespace Eigen; int main() { MatrixXd m = MatrixXd::Random(3,3); // 初始化动态数组 m = (m + MtrixXd::Constant(3,3,1.2)) * 50; cout<<"m = "<<endl<<m<<endl; VectorXd v(3); v<<1,2,3; // 逗号初始化 cout<< "m * v = "<< endl << m * v << endl; }
2.3 comma-initializer
RowVectorXd vec1(3); vec1 << 1,2,3; std::cout<<"vec1 = "<< vec1 << std::endl; RowVectorXd vec2(4); vec1 << 1,2,3,4; std::cout<<"vec2 = "<< vec2 << std::endl; RowVectorXd joined(7); joined << vec1,vec2; std::cout<<"joined = "<< joined << std::endl;
2.4 常用的初始化方法(初始化为0,初始化为1,初始化为单位矩阵)
#include<iostream> #include<Eigen/Dense> using namespace std; using namespace Eigen; int main() { MatrixXd m0 = MatrixXd::Random(3,3); // 随机初始化,值在(-1,1) MatrixXd m1 = MatrixXd::Constant(3,3,2.4); // 常量初始化,值全部是2.4 MatrixXd m2 = MatrixXd::Zero(); // 零初始化 MatrixXd m3 = MatrixXd::Ones(); // 初始化1 MatrixXd m4 = MatrixXd:Identity();// 初始化为单位矩阵 }
2.5 调整矩阵大小
矩阵的当前大小可以通过rows()、cols()和size()获取。调整动态大小矩阵的大小是通过resize()方法实现的。 动态矩阵可以随意调整矩阵大小,固定尺寸的矩阵无法调整大小。
MatrixXd m0(x,x); m0.resize(4,4); VectorXd v(3); v.resize(5);
2.6 固定尺寸与动态尺寸的选择
尽可能使用固定尺寸来显示非常小的尺寸,尽可能使用动态尺寸显示较大的尺寸。对于小尺寸,尤其是对于小于(大约)16的尺寸,使用固定尺寸对性能有极大的好处。因为它使Eigen避免了动态内存分配并展开循环。在内部,固定大小的本征矩阵只是一个金蛋的数组.
3 矩阵类
在Eigen中,所有matrices和vectors都是Matrix模板类的对象。vectors只是matrices的一种特殊情况,具有1行和1列。
3.1 Matrix的前三个模板参数
该矩阵类需要留个模板参数,但一般只需了解前三个参数即可。剩下的三个参数具有默认值。 Matrix的三个必需模板参数: Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> Scalar是标量类型,即系数的类型。 RowsAtCompileTime/ColsAtCompileTime是在编译时已知矩阵的行数和列数。
typedef Matrix<int, 1, 2> RowVector2i; typedef Matrix<double,Dynamic,Dynamic> MatrixXd; typedef Matrix<int, Dynamic, 1> VectorXi; typedef Array<float, Dynamic, Dynamic> ArrayXXf; typedef Array<double, Dynamic,1> ArrayXd; typedef Array<int, 1, Dynamic> RowArrayXi; typedef Array<float, 3, 3> Array33f; typedef Array<float, 4,1> Array4f; ...
3.2 vectors
在Eigen中,vectors只是Matrix的一种特殊情况,具有1行或者1列。他们只有1列的情况最为常见,这样的向量称为列向量,通常缩写为向量。只有1行的向量称为行向量。
3.3 动态矩阵
动态矩阵在编译的时候不知道其大小,需要在运行的时候才能确定其大小。 typedef Matrix<dobule, Dynamic, Dynamic> MatrixXd;
4 Array类的介绍
Eigen不仅提供了Matrix和Vector结构,还提供了Array结构。区别如下:Matrix和Vector就是线性代数中定义的矩阵和向量,所有的数学运算都和数学上一致。但是存在一个问题就是数学上的定义不一定完全满足实际要求。比如,数学上并没有定义一个矩阵和一个标量的加法运算。但是如果我们想给一个矩阵的每个元素都加上同一个数,那么这个操作就需要我们自己去实现,这显然并不方便。 Array提供了一个Array类,为我们提供了大量的矩阵未定义的操作,且Array和Matrix之间很容易相互转换,所以相当于给矩阵提供更多的方法,也为使用者的不同需求提供了更多的选择。 Array<typename Scalar, int row, int cols>
4.1 Array初始化,加减乘除操作
Eigen::Array类重载了+-*/运算符,可以直接用这些运算符对Array对象进行操作,相乘操作是对应的数字相乘,相除操作是对应的元素相除。
ArrayXXf a(3,3); ArrayXXf b(,3,3); a<<1,2,3,4,5,6,7,8,9; b<<1,2,3,1,2,3,1,2,3; cout<<"a + b = "<<endl<<a + b<<endl; cout<<"a - 2 = "<<endl<<a - 2<<endl; // 所有元素减2 cout<<"a * b = "<<endl <<a * b<<endl; // 对应元素相乘 cout<<"a / b = "<<endl <<a / b<<endl; // 对应元素相除
4.2 Array类的其它操作
Eigen::Array类还定义了绝对值abs(),开平方根sqrt()及找对应元素最小值操作min()。
ArrayXXf a = ArrayXXf::Random(3,3); a *= 2; cout<<"a = "<<endl<<a<<endl; cout<<"a.abs() = "<<endl<<a.abs()<<endl; cout<<"a.abs().sqrt()"<<endl<<a.abs().sqrt()<<endl; cout<<"a.min(a.abs().sqrt()) = "<<endl<<a.min(a.abs().sqrt())<<endl;
6. Matrix和Array之间的相互转换
Matrix类和Array类之间可以相互转换,必须显示转换,才能对他们进行加减乘除运算。
Array44f a1,a2; Matrix4f m1,m2; m1 = a1 * a2; a1 = m1 * m2; a2 = a1 + m1.array(); m2 = a1.matrix() + m1; ArrayWrapper<Matrix4f> m1a(m1);// m1a是m1.array()的别名,他们共享相同的系数 MatrixWrapper<Array44f>a1m(a1);
7. 矩阵转置、共轭、共轭转置
下面介绍矩阵的一些操作。
7.1 转置与共轭
对矩阵的转置、共轭和共轭转置由成员函数transpose()/conjugate()/adjoint()实现。
MatrixXcf a = MatrixXcf::Random(2,2); cout<<"a = "<<endl <<a <<endl; cout<<"a.transpose = "<<endl<<a.transpose()<<endl; cout<<"a.conjugate = "<<endl<<a.conjugate()<<endl; cout<<"a.adjoint = "<<endl<<a.adjoint()<<endl;
7.2 转置需要注意的事项
a = a.tanspose(); 无法运行,这称为别名问题。在Debug模式下当assertions没有禁止时,这种问题会被自动检测到。要避免错误,可以使用in-place转置。
Matrix2i a; a<<1,2,3,4; cout<<" a = \n"<<a<<endl; // a = a.transpose();// 不要这样写代码,无法运行 a.transposeInPlace(); cout<<"a.transpose = \n"<<a<<endl;
8. 点积和叉积
对于点积和叉积,直接使用dot()和cross()方法。
Vector3d v(1,2,3); Vector3d w(0,1,2); cout<<"Dot product: "<<v.dot(w)<<endl; double dp = v.adjoint() * w; cout<<"Dot product via a matrix product:"<<dp <<endl; cout<<"Cross product:\n"<<v.cross(w) <<endl;
注意:叉积仅仅用于尺寸为3的向量!点积可以用于任意尺寸的向量,当时用复数时,Eigen的点积作为第一个变量为共轭线性的,第二个为线性的。
9. 矩阵的基础的算术(求和、平均值等)
Eigen提供了一些对于矩阵或向量的规约操作,如sum()、prod()、maxCoeff()、minCoeff()
Eigen::Matrix2d mat; mat<<1,2,3,4; cout<<"mat.sum: "<<mat.sum() <<endl; cout<<"mat.prod: "<<mat.prod() <<endl; cout<<"mat.mean: "<<mat.mean() <<endl; cout<<"mat.minCoeff: "<<mat.minCoeff() <<endl; cout<<"mat.maxCoeff: "<<mat.maxCoeff() <<endl; cout<<"mat.trace: "<<mat.trace() <<endl;
10. Eigen块操作
10.1 块基本操作
块指的是矩阵或者数组中的一个矩形区域,块表达式可以用于左值或者右值,同样不会耗费运行时间,由编译器优化。 Eigen中最常用的快操作是block()方法,共有两个版本
索引从0开始,两个版本都可用于固定尺寸或者动态尺寸的矩阵和数组。这两个表达式语义上相同,唯一的区别是如果块的尺寸比较小的话固定尺寸版本的块操作运行更快,但是需要在便器阶段知道大小。
#include <Eigen/Dense> #include <iostream> using namespace std; int main() { Eigen::MatrixXf m(4, 4); // 初始化m矩阵 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { m(i, j) = j + 1 + i * 4; } } cout << "m: " << endl << m << endl; cout << "Block in the middle" << endl; cout << m.block<2, 2>(1, 1) << endl << endl;// m.block<i,j> (a,b) 表示从第(a+1)行(b+1)列开始,截图i行,j列 for (int i = 1; i <= 3; ++i) { cout << "Block of size " << i << "x" << i << endl; cout << m.block(0, 0, i, i) << endl << endl;// m.block(a,b,i,j) 表示从第(a+1)行(b+1)列开始,截图i行,j列 } }
注意:m.block<i,j>(a,b)表示从第(a+1)行(b+1)列开始,截i行j列。 上述例子中的块操作方法作为表达式的右值,意味着是只读形式的,然而,快操作也可以作为左值使用,意味着可以给其幅值,下面的例子说明了这一点,当然对于矩阵的操作是一样的。
#include <Eigen/Dense> #include <iostream> using namespace std; using namespace Eigen; int main() { Array22f m; m << 1, 2, 3, 4; Array44f a = Array44f::Constant(0.6); cout << "Here is the array a:" << endl << a << endl << endl; a.block<2, 2>(1, 1) = m; cout << "Here is now a with m copied into its central 2x2 block:" << endl << a << endl << endl; a.block(0, 0, 2, 3) = a.block(2, 1, 2, 3); cout << "Here is now a with bottom-right 2x3 block copied into top-left 2x2 block:" << endl << a << endl << endl; }
尽管block()方法可用于任何形式的块操作,还是有特殊的方法用于特殊情况的,主要是为了更好的性能。说到性能,最重要的是在编译阶段给EIgen尽可能多的信息。比如,如果你的块是一个矩阵中的一列,那么使用col()方法会更好。
10.2 行和列
行和列是一种特殊的块,Eigen提供了特殊的方法:col()和row()。
#include <Eigen/Dense> #include <iostream> using namespace std; int main() { Eigen::MatrixXf m(4, 4); // 数组初始化 for (int i = 0; i < 4; ++i) { for (int j = 0; j < 4; ++j) { m(i, j) = j + 1 + i * 4; } } 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; }
10.3 边角相关的操作
Eigen同样提供了对于挨着矩阵或数组的边、角的特殊操作方法,比如topLeftCorner()方法可用于操作矩阵左上角的区域。总结如下:
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<2>() =" << endl << m.bottomRows<2>() << endl << endl; m.topLeftCorner(1, 3) = m.bottomRightCorner(3, 1).transpose(); cout << "After assignment, m = " << endl << m << endl;
10.4 对于向量的块操作
Eigen也提供了一些针对向量和一维数组的块操作方法:
Eigen::ArrayXf v(6); v << 1, 2, 3, 4, 5, 6; cout << "v.head(3) =" << endl << v.head(3) << endl << endl; cout << "v.tail<3>() = " << endl << v.tail<3>() << endl << endl; v.segment(1,4) *= 2; cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
11 范数计算
向量的平方范数由squaredNorm()获得,等价于向量对自身做点积,也等同于所有元素额平方和。Eigen也提供了norm()范数,返回的是squaredNorm()的根。这些操作也适用于矩阵。如果想使用其他元素级的范数,使用lpNorm()方法,当求无穷范数时,模板参数p可以取特殊值Infinity,得到的是所有元素的最大绝对值。
#include <Eigen/Dense> #include <iostream> using namespace std; using namespace Eigen; int main() { VectorXf v(2); MatrixXf m(2, 2), n(2, 2); v << -1, 2; m << 1, -2, -3, 4; cout << "v.squaredNorm() = " << v.squaredNorm() << endl; cout << "v.norm() = " << v.norm() << endl; cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl; cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl; cout << endl; cout << "m.squaredNorm() = " << m.squaredNorm() << endl; cout << "m.norm() = " << m.norm() << endl; cout << "m.lpNorm<1>() = " << m.lpNorm<1>() << endl; cout << "m.lpNorm<Infinity>() = " << m.lpNorm<Infinity>() << endl; }
矩阵的1范数和无穷范数也可以用下面的方法计算:
#include <Eigen/Dense> #include <iostream> using namespace Eigen; using namespace std; int main() { MatrixXf m(2,2); m << 1,-2, -3,4; cout << "1-norm(m) = " << m.cwiseAbs().colwise().sum().maxCoeff() << " == " << m.colwise().lpNorm<1>().maxCoeff() << endl; cout << "infty-norm(m) = " << m.cwiseAbs().rowwise().sum().maxCoeff() << " == " << m.rowwise().lpNorm<1>().maxCoeff() << endl; }
12 布尔规约
如下的操作得到的是布尔值 all()返回真,如果矩阵或数组的所有元素为真 any()返回真,如果矩阵或数组至少有一个元素为真 count()返回元素为真的个数
#include <Eigen/Dense> #include <iostream> using namespace std; using namespace Eigen; int main() { ArrayXXf a(2,2); a << 1,2,3,4; cout << "(a > 0).all() = " << (a > 0).all() << endl; cout << "(a > 0).any() = " << (a > 0).any() << endl; cout << "(a > 0).count() = " << (a > 0).count() << endl; cout << endl; cout << "(a > 2).all() = " << (a > 2).all() << endl; cout << "(a > 2).any() = " << (a > 2).any() << endl; cout << "(a > 2).count() = " << (a > 2).count() << endl; }
13迭代
当需要获得元素在矩阵或数组中的位置时使用迭代。
#include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { Eigen::MatrixXf m(2,2); m << 1, 2,3, 4; //get location of maximum MatrixXf::Index maxRow, maxCol; float max = m.maxCoeff(&maxRow, &maxCol); //get location of minimum MatrixXf::Index minRow, minCol; float min = m.minCoeff(&minRow, &minCol); cout << "Max: " << max << ", at: " << maxRow << "," << maxCol << endl; cout << "Min: " << min << ", at: " << minRow << "," << minCol << endl; }
14 部分规约
部分规约指的是对矩阵或数组按行或列进行的操作,对每一列或者行进行规约操作时得到的是一个列或者行向量。如下例子得到矩阵每一列的最大值并存入一个行向量中。
#include <iostream> #include <Eigen/Dense> using namespace std; int main() { Eigen::MatrixXf mat(2,4); mat << 1, 2, 6, 9, 3, 1, 7, 2; std::cout << "Column's maximum: \n" <<mat.colwise().maxCoeff() << std::endl; }
同样也可以得到每一行的最大值,返回一个列向量。
#include <iostream> #include <Eigen/Dense> 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; }
15 广播机制
广播的概念类似于部分规约,不同之处在于广播通过对向量在一个方向上的复制,将向量解释成矩阵。如下例子将一个列向量加到矩阵的每一列中。
#include <iostream> #include <Eigen/Dense> 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; std::cout << "Broadcasting result: " << std::endl; std::cout << mat << std::endl; }
可以将mat.colwise()+=v理解成两种等价的方式,它将列向量加到矩阵的每一列中;或者将列向量复制4次的得到一个2x4的矩阵,之后进行矩阵的相加运算:
+=、+和-运算符也可以按列或行操作。在数组中也可以用=、/=、和/运算符执行元素级的按行或列乘除运算。但不能用在矩阵上,如果想用v(0)乘以矩阵的第0列,v(1)乘以矩阵的第1列…使用mat = matv.asDiagonal()。
结合广播和其他操作
广播也可以和其他操作结合,比如矩阵或数组的运算、规约和部分规约操作。下面介绍一个更加复杂的例子,演示了在矩阵中找到和给定向量最接近的一列,使用到了欧氏距离。
#include <iostream> #include <Eigen/Dense> using namespace std; using namespace Eigen; int main() { Eigen::MatrixXf m(2,4); Eigen::VectorXf v(2); m << 1, 23, 6, 9, 3, 11, 7, 2; v << 2, 3; MatrixXf::Index index; // find nearest neighbour (m.colwise() - v).colwise().squaredNorm().minCoeff(&index); cout << "Nearest neighbour is column " << index << ":" << endl; cout << m.col(index) << endl; }