入门
安装
下载并解压Eigen源代码,并把里面的Eigen子目录放到包含目录里即可
第一个简单的程序
#include <iostream>
#include <Eigen/Dense>
using Eigen::MatrixXd;
int main()
{
MatrixXd m(2,2);
m(0,0) = 3;
m(1,0) = 2.5;
m(0,1) = -1;
m(1,1) = m(1,0) + m(0,1);
std::cout << m << std::endl;
}
例二:矩阵和矢量
运行时设定矩阵/向量大小
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
MatrixXd m = MatrixXd::Random(3,3);
// 当元素类型为整型时,范围为整型的范围;当元素类型为浮点型时,范围为[-1,1]
m = (m + MatrixXd::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;
}
编译时设定矩阵/向量大小
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
Matrix3d m = Matrix3d::Random();
m = (m + Matrix3d::Constant(1.2)) * 50;
cout << "m =" << endl << m << endl;
Vector3d v(1,2,3);
cout << "m * v =" << endl << m * v << endl;
}
两种方式的选择原则:4*4或更小的矩阵用后者,否则用前者
矩阵类
矩阵模板的前三个参数(必填)
Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
- Scalar:矩阵元素的类型(float、double等)
- RowsAtCompileTime、ColsAtCompileTime:编译时的行数、列数
向量
typedef Matrix<float, 3, 1> Vector3f;
typedef Matrix<int, 1, 2> RowVector2i;
动态大小
typedef Matrix<double, Dynamic, Dynamic> MatrixXd;
typedef Matrix<int, Dynamic, 1> VectorXi;
构造函数
3*3,元素未初始化
Matrix3f a;
动态大小,0*0,暂无元素
MatrixXf b;
动态大小,10*15,元素未初始化
MatrixXf a(10,15);
动态大小,30*1,元素未初始化
VectorXf b(30);
固定大小的构造函数可以接受大小,但不会起作用
Matrix3f a(3,3);
初始化元素
Vector2i a(1, 2); // A column vector containing the elements {1, 2}
Matrix<int, 5, 1> b {1, 2, 3, 4, 5}; // A row-vector containing the elements {1, 2, 3, 4, 5}
Matrix<int, 1, 5> c = {1, 2, 3, 4, 5}; // A column vector containing the elements {1, 2, 3, 4, 5}
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},
};
VectorXd a {{1.5, 2.5, 3.5}}; // A column-vector with 3 coefficients
RowVectorXd b {{1.0, 2.0, 3.0, 4.0}}; // A row-vector with 4 coefficients
访问元素
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
MatrixXd m(2, 2);
m(0, 0) = 3;
m(1, 0) = 2.5;
m(0, 1) = -1;
m(1, 1) = m(1, 0) + m(0, 1);
std::cout << "Here is the matrix m:\n" << m << std::endl;
std::cout << m(2) << std::endl; // column-major
VectorXd v(2);
v(0) = 4;
v(1) = v(0) - 1;
std::cout << "Here is the vector v:\n" << v << std::endl;
}
逗号初始化
Matrix3f m;
m << 1, 2, 3,
4, 5, 6,
7, 8, 9;
std::cout << m;
改变大小
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
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;
}
resize()不能改变确定大小的矩阵,但其他函数可以正常工作
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix4d m;
m.resize(4,4); // no operation
std::cout << "The matrix m is of size "
<< m.rows() << "x" << m.cols() << std::endl;
}
矩阵赋值
左侧矩阵会变成右侧矩阵的大小
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;
固定大小与动态大小
如何选择
- 固定大小:小矩阵(小于16),能用的时候
- 动态大小:大矩阵(大于16),不得不用的时候
选填的模板参数
Matrix<typename Scalar,
int RowsAtCompileTime,
int ColsAtCompileTime,
int Options = 0,
int MaxRowsAtCompileTime = RowsAtCompileTime,
int MaxColsAtCompileTime = ColsAtCompileTime>
- Options:决定存储方式,默认为0(ColMajor),可填1(RawMajor)
- MaxRowsAtCompileTime、MaxColsAtCompileTime:定义最大行数、列数,此时即使前面定义为Dynamic,也可以避免动态内存分配
快捷类型定义
- MatrixNt:Matrix<type, N, N>
- VectorNt:Matrix<type, N, 1>
- RowVectorNt:Matrix<type, 1, N>
其中
- N可以是2、3、4、X
- t可以是i(整型)、f(浮点型)、d(双精度浮点型)、cf(浮点型复数)、cd(双精度浮点型复数)
矩阵和向量运算
加减
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
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;
Vector3d v(1, 2, 3);
Vector3d w(1, 0, 0);
std::cout << "-v + w - v =\n" << -v + w - v << std::endl;
}
标量乘除
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d a;
a << 1, 2,
3, 4;
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;
}
转置和共轭
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;
transpose()和adjoint()的一点问题
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;
矩阵相乘
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{
Matrix2d mat;
mat << 1, 2,
3, 4;
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;
}
点积和叉积
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
Vector3d v(1,2,3);
Vector3d w(0,1,2);
cout << "Dot product: " << v.dot(w) << endl;
double dp = v.adjoint()*w; // automatic conversion of the inner product to a scalar
cout << "Dot product via a matrix product: " << dp << endl;
cout << "Cross product:\n" << v.cross(w) << endl;
}
叉积只能用于大小为3的矢量,点积可用于大小任意的矢量
简便操作
#include <iostream>
#include <Eigen/Dense>
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;
操作可行性检查
Matrix3f m;
Vector4f v;
v = m*v; // Compile-time error: YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES
MatrixXf m(3,3);
VectorXf v(4);
v = m * v; // Run-time assertion failure here: "invalid matrix product"
数组类
数组模板
声明与矩阵一致
Array<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime>
Type | Typedef |
Array<float,Dynamic,1> | ArrayXf |
Array<float,3,1> | Array3f |
Array<double,Dynamic,Dynamic> | ArrayXXd |
Array<double,3,3> | Array33d |
访问数组元素
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
ArrayXXf m(2,2);
// assign some values coefficient by coefficient
m(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 output
cout << m << endl << endl;
// using the comma-initializer is also allowed
m << 1.0,2.0,
3.0,4.0;
// print values to standard output
cout << m << endl;
}
加减
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
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;
// Adding two arrays
cout << "a + b = " << endl << a + b << endl << endl;
// Subtracting a scalar from an array
cout << "a - 2 = " << endl << a - 2 << endl;
}
乘法
对应位置元素相乘
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
ArrayXXf a(2,2);
ArrayXXf b(2,2);
a << 1,2,
3,4;
b << 5,6,
7,8;
cout << "a * b = " << endl << a * b << endl;
}
其他元素操作
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
int main()
{
ArrayXf a = ArrayXf::Random(5);
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;
}
矩阵与数组之间转换
#include <Eigen/Dense>
#include <iostream>
using namespace Eigen;
using namespace std;
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;
cout << "-- Matrix m*n: --" << endl << result << endl << endl;
result = m.array() * n.array();
cout << "-- Array m*n: --" << endl << result << endl << endl;
result = m.cwiseProduct(n);
cout << "-- With cwiseProduct: --" << endl << result << endl << endl;
result = m.array() + 4;
cout << "-- Array m + 4: --" << endl << result << endl << endl;
}
分块操作
使用分块操作
操作 | 动态大小 | 固定大小 |
从i、j开始,大小为p、q的分块 | matrix.block(i,j,p,q); | matrix.block<p,q>(i,j); |
当分块较小时,使用固定大小代码执行更快,但是需要在编译时知道分块的大小(参数为数字而不能是变量)
#include <Eigen/Dense>
#include <iostream>
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 <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 2x3 block:" << endl << a << endl << endl;
}
行与列
#include <Eigen/Dense>
#include <iostream>
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;
}
边角处的分块
动态大小 | 固定大小 |
matrix.topLeftCorner(p,q); | matrix.topLeftCorner<p,q>(); |
matrix.bottomLeftCorner(p,q); | matrix.bottomLeftCorner<p,q>(); |
matrix.topRightCorner(p,q); | matrix.topRightCorner<p,q>(); |
matrix.bottomRightCorner(p,q); | matrix.bottomRightCorner<p,q>(); |
matrix.topRows(q); | matrix.topRows<q>(); |
matrix.bottomRows(q); | matrix.bottomRows<q>(); |
matrix.leftCols(p); | matrix.leftCols<p>(); |
matrix.rightCols(q); | matrix.rightCols<q>(); |
matrix.middleCols(i,q); | matrix.middleCols<q>(i); |
matrix.middleRows(i,q); | matrix.middleRows<q>(i); |
#include <Eigen/Dense>
#include <iostream>
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<2>() =" << endl << m.bottomRows<2>() << endl << endl;
m.topLeftCorner(1,3) = m.bottomRightCorner(3,1).transpose();
cout << "After assignment, m = " << endl << m << endl;
}
矢量的分块操作
动态大小 | 固定大小 |
vector.head(n); | vector.head<n>(); |
vector.tail(n); | vector.tail<n>(); |
vector.segment(i,n); | vector.segment<n>(i); |
#include <Eigen/Dense>
#include <iostream>
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<3>() = " << endl << v.tail<3>() << endl << endl;
v.segment(1,4) *= 2;
cout << "after 'v.segment(1,4) *= 2', v =" << endl << v << endl;
}
切片和索引
概要
DenseBase::operator()(const RowIndices&, const ColIndices&)
参数可以是
- 整数(变量)
- Eigen::all、Eigen::last
- Eigen::seq、Eigen::seqN、Eigen::lastN
- 1为向量或数组(Eigen、std、C数组都可以)
基础切片
seq(firstIdx,lastIdx) | seq(2,5) <=> {2,3,4,5} |
seq(firstIdx,lastIdx,incr) | seq(2,8,2) <=> {2,4,6,8} |
seqN(firstIdx,size) | seqN(2,5) <=> {2,3,4,5,6} |
seqN(firstIdx,size,incr) | seqN(2,3,3) <=> {2,5,8} |
firstIdx和lastIdx可以是Eigen::last
切片 | 等效的分块 |
v(lastN(n)) | v.tail(n) |
v(lastN(m), lastN(n)) | A.bottomRightCorner(m,n) |
v(lastN(m), lastN(n)) | A.bottomRightCorner(m,n) |
A(all, lastN(n,3)) |
提高性能
v(seq(last-fix<7>, last-fix<2>))
v(seqN(last-7, fix<6>))
A(all, seq(0,last,fix<2>))
倒序切片
A(all, seq(20, 10, fix<-2>))
A(seqN(last, n, fix<-1>), all)
A(lastN(n).reverse(), all)
索引数组
std::vector<int> ind{4,2,5,5,3};
MatrixXi A = MatrixXi::Random(4,6);
cout << "Initial matrix A:\n" << A << "\n\n";
cout << "A(all,ind):\n" << A(all,ind) << "\n\n";
#if EIGEN_HAS_STATIC_ARRAY_TEMPLATE
MatrixXi A = MatrixXi::Random(4,6);
cout << "Initial matrix A:\n" << A << "\n\n";
cout << "A(all,{4,2,5,5,3}):\n" << A(all,{4,2,5,5,3}) << "\n\n";
#endif
ArrayXi ind(5); ind<<4,2,5,5,3;
MatrixXi A = MatrixXi::Random(4,6);
cout << "Initial matrix A:\n" << A << "\n\n";
cout << "A(all,ind-1):\n" << A(all,ind-1) << "\n\n";
自定义索引列表(好像懂了)
struct pad {
Index size() const { return out_size; }
Index operator[] (Index i) const { return std::max<Index>(0,i-(out_size-in_size)); }
Index in_size, out_size;
};
Matrix3i A;
A.reshaped() = VectorXi::LinSpaced(9,1,9);
cout << "Initial matrix A:\n" << A << "\n\n";
MatrixXi B(5,5);
B = A(pad{3,5}, pad{3,5});
cout << "A(pad{3,N}, pad{3,N}):\n" << B << "\n\n";
高级初始化
逗号初始化
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;
特殊矩阵和数组
零
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";
常量
Matrix3d::Constant(value)
MatrixXd::Constant(rows, cols, value)
随机数
Matrix3d::Random()
MatrixXd::Random(rows, cols)
单位矩阵(仅适用于矩阵)
Matrix3d::Identity()
MatrixXd::Identity(rows, cols)
元素线性分布
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;
简便操作
- setZero()
- setIdentity()
- setLinSpaced()
用作临时变量
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;
int main()
{
MatrixXd m = MatrixXd::Random(3,3);
m = (m + MatrixXd::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;
}
MatrixXf mat = MatrixXf::Random(2, 3);
std::cout << mat << std::endl << std::endl;
mat = (MatrixXf(2,2) << 0, 1, 1, 0).finished() * mat;
std::cout << mat << std::endl;
注意finish()不能忘
Reductions, visitors and broadcasting(不会翻译~)
Reductions
#include <iostream>
#include <Eigen/Dense>
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;
}
计算范数
#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;
}
布尔操作
#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;
}
Visitors
Eigen::index <=> std::ptrdiff_t
#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;
}
局部reductions
#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: " << std::endl
<< 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;
}
将局部reductions与其他操作混合
#include <iostream>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;
int main()
{
MatrixXf mat(2,4);
mat << 1, 2, 6, 9,
3, 1, 7, 2;
MatrixXf::Index maxIndex;
float maxNorm = mat.colwise().sum().maxCoeff(&maxIndex);
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;
}
广播
#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;
}
#include <iostream>
#include <Eigen/Dense>
using namespace std;
int main()
{
Eigen::MatrixXf mat(2,4);
Eigen::VectorXf v(4);
mat << 1, 2, 6, 9,
3, 1, 7, 2;
v << 0,1,2,3;
//add v to each row of m
mat.rowwise() += v.transpose();
std::cout << "Broadcasting result: " << std::endl;
std::cout << mat << std::endl;
}
把broadcasting与其他操作混合
#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;
}
Reshape
Reshaped 2D views
Matrix4i m = Matrix4i::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is m.reshaped(2, 8):" << endl << m.reshaped(2, 8) << endl;
按column-major把原矩阵的元素一个个摆进新矩阵
1D linear views
Matrix4i m = Matrix4i::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is m.reshaped().transpose():" << endl << m.reshaped().transpose() << endl;
cout << "Here is m.reshaped<RowMajor>().transpose(): " << endl << m.reshaped<RowMajor>().transpose() << endl;
ReshapeInPlace
MatrixXi m = Matrix4i::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is m.reshaped(2, 8):" << endl << m.reshaped(2, 8) << endl;
m.resize(2,8);
cout << "Here is the matrix m after m.resize(2,8):" << endl << m << endl;
Resize的ColMajor/RawMajor取决于矩阵的存储方式
Matrix<int,Dynamic,Dynamic,RowMajor> m = Matrix4i::Random();
cout << "Here is the matrix m:" << endl << m << endl;
cout << "Here is m.reshaped(2, 8):" << endl << m.reshaped(2, 8) << endl;
cout << "Here is m.reshaped<AutoOrder>(2, 8):" << endl << m.reshaped<AutoOrder>(2, 8) << endl;
m.resize(2,8);
cout << "Here is the matrix m after m.resize(2,8):" << endl << m << endl;
由于aliasing
A = A.reshaped(2,8); // 不允许
A = A.reshaped(2,8).eval(); // 允许
STL迭代器和算法
遍历一维数组和向量
VectorXi v = VectorXi::Random(4);
cout << "Here is the vector v:\n";
for(auto x : v) cout << x << " ";
cout << "\n";
Array4i v = Array4i::Random().abs();
cout << "Here is the initial vector v:\n" << v.transpose() << "\n";
std::sort(v.begin(), v.end());
cout << "Here is the sorted vector v:\n" << v.transpose() << "\n";
遍历二维数组和矩阵的元素
Matrix2i A = Matrix2i::Random();
cout << "Here are the coeffs of the 2x2 matrix A:\n";
for(auto x : A.reshaped())
cout << x << " ";
cout << "\n";
遍历二维数组和矩阵的行或者列
ArrayXXi A = ArrayXXi::Random(4,4).abs();
cout << "Here is the initial matrix A:\n" << A << "\n";
for(auto row : A.rowwise())
std::sort(row.begin(), row.end());
cout << "Here is the sorted matrix A:\n" << A << "\n";
与C/C++数组进行交互:Map类
声明
Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >
例子
Map<MatrixXf> mf(pf,rows,columns); // pf是float数组的指针
Map<const Vector4i> mi(pi); // pi是int数组的指针
高级用法(没懂)
Map<typename MatrixType, int MapOptions, typename StrideType>
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;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;
使用Map变量(有点意义不明?)
typedef Matrix<float,1,Dynamic> MatrixType;
typedef Map<MatrixType> MapType;
typedef Map<const MatrixType> MapTypeConst; // a read-only map
const 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 m2
MapType m2map(p,m2.size()); // m2map shares data with m2
MapTypeConst m2mapconst(p,m2.size()); // a read-only accessor for m2
cout << "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 array
cout << "Updated m2: " << m2 << endl;
cout << "m2 coefficient 2, constant accessor: " << m2mapconst(2) << endl;
/* m2mapconst(2) = 5; */ // this yields a compile-time error
更改映射数组
int data[] = {1,2,3,4,5,6,7,8,9};
Map<RowVectorXi> v(data,4);
cout << "The mapped vector v is: " << v << "\n";
new (&v) Map<RowVectorXi>(data+4,5);
cout << "Now v is: " << v << "\n";
后面这段完全没懂
Map<Matrix3f> A(NULL); // don't try to use this matrix yet!
VectorXf b(n_matrices);
for (int i = 0; i < n_matrices; i++)
{
new (&A) Map<Matrix3f>(get_matrix_pointer(i));
b(i) = A.trace();
}
Aliasing
赋值语句左右两边有同一个变量,有可能导致奇怪问题
例子
MatrixXi mat(3,3);
mat << 1, 2, 3, 4, 5, 6, 7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
// This assignment shows the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);
cout << "After the assignment, mat = \n" << mat << endl;
解决方法:eval()
MatrixXi mat(3,3);
mat << 1, 2, 3, 4, 5, 6, 7, 8, 9;
cout << "Here is the matrix mat:\n" << mat << endl;
// The eval() solves the aliasing problem
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2).eval();
cout << "After the assignment, mat = \n" << mat << endl;
xxxInPlace()比xxx().eval()更好,有的话尽量用
Aliasing和矩阵/数组运算
左值的元素仅取决于右值的对应元素时,不会出现aliasing
MatrixXf mat(2,2);
mat << 1, 2, 4, 7;
cout << "Here is the matrix mat:\n" << mat << endl << endl;
mat = 2 * mat;
cout << "After 'mat = 2 * mat', mat = \n" << mat << endl << endl;
mat = mat - MatrixXf::Identity(2,2);
cout << "After the subtraction, it becomes\n" << mat << endl << endl;
ArrayXXf arr = mat;
arr = arr.square();
cout << "After squaring, it becomes\n" << arr << endl << endl;
// Combining all operations in one statement:
mat << 1, 2, 4, 7;
mat = (2 * mat - MatrixXf::Identity(2,2)).array().square();
cout << "Doing everything at once yields\n" << mat << endl << endl;
Aliasing和矩阵乘法
矩阵乘法默认会进行eval(),如果左值没有在右值中出现,可以用noalias()加速
MatrixXf matA(2,2), matB(2,2);
matA << 2, 0, 0, 2;
// Simple but not quite as efficient
matB = matA * matA;
cout << matB << endl << endl;
// More complicated but also more efficient
matB.noalias() = matA * matA;
cout << matB;
例外:矩阵的形状发生变化,且矩阵乘积没有之间赋值给左值,此时需要手动eval()
总结
当赋值语句左边和右边出现同一个矩阵或者数组,则有可能发生aliasing
- 如果左边的元素仅取决于右边的对应元素,则aliasing不影响
- 矩阵乘法默认会用eval(),可以用noalias()加速
- 在其他情况下,需要用eval(),xxxInPlace()来避免aliasing
存储方式
column-major和row-major
描述矩阵的元素是按什么顺序存储的
- column-major:从左到右、从上到下(默认)
- row-major:从上到下、从左到右
Matrix<int, 3, 4, ColMajor> Acolmajor;
Acolmajor << 8, 2, 2, 9,
9, 1, 4, 4,
3, 5, 4, 5;
cout << "The matrix A:" << endl;
cout << Acolmajor << endl << endl;
cout << "In memory (column-major):" << endl;
for (int i = 0; i < Acolmajor.size(); i++)
cout << *(Acolmajor.data() + i) << " ";
cout << endl << endl;
Matrix<int, 3, 4, RowMajor> Arowmajor = Acolmajor;
cout << "In memory (row-major):" << endl;
for (int i = 0; i < Arowmajor.size(); i++)
cout << *(Arowmajor.data() + i) << " ";
cout << endl;
两种方式的矩阵或数组之间可以相互赋值、运算
如何选择
- 根据应用需求、其他库的需求
- 遍历行的算法用row-major、遍历列的算法用column-major理论上更快
- 没有特殊需求就用默认的column-major
关于未对齐数组的assertion的解释
程序运行过程中可能出现下面的报错
my_program: path/to/eigen/Eigen/src/Core/DenseStorage.h:44:
Eigen::internal::matrix_array<T, Size, MatrixOptions, Align>::internal::matrix_array()
[with T = double, int Size = 2, int MatrixOptions = 2, bool Align = true]:
Assertion `(reinterpret_cast<size_t>(array) & (sizemask)) == 0 && "this assertion
is explained here: http://eigen.tuxfamily.org/dox-devel/group__TopicUnalignedArrayAssert.html
READ THIS WEB PAGE !!! ****"' failed.
使用新的编译器可以解决这个问题,如GCC>=7, clang>=5, MSVC>=19.12
该报错不能反映程序那里出现问题,需要手动逐步运行找
四个导致该问题的原因
仅针对固定大小的Eigen变量
- Eigen变量作为结构体的成员
- STL容器或手动内存分配
- 使用Eigen变量传值
- 编译器对堆栈对齐做出错误假设(例如Windows上的GCC)
解决问题的可能方法
- 对出问题的变量使用DontAlign
- 将EIGEN_MAX_STATIC_ALIGN_BYTES定义为0
- 定义EIGEN_DONT_VECTORIZE和EIGEN_DISABLE_UNALIGNED_ARRAY_ASSERT
固定大小的可矢量化的Eigen变量
解释:有固定大小、且大小是16字节的整数倍的变量
- Eigen::Vector2d
- Eigen::Vector4d
- Eigen::Vector4f
- Eigen::Matrix2d
- Eigen::Matrix2f
- Eigen::Matrix4d
- Eigen::Matrix4f
- Eigen::Affine3d
- Eigen::Affine3f
- Eigen::Quaterniond
- Eigen::Quaternionf
Eigen变量作为结构体成员
class Foo
{
...
Eigen::Vector4d v;
...
public:
EIGEN_MAKE_ALIGNED_OPERATOR_NEW
};
...
Foo *foo = new Foo;
使用STL容器
把
std::map<int, Eigen::Vector4d>
改成
std::map<int, Eigen::Vector4d, std::less<int>, Eigen::aligned_allocator<std::pair<const int, Eigen::Vector4d> > >
把Eigen变量传值到函数里面
把
void my_function(Eigen::Vector2d v);
改成
void my_function(const Eigen::Vector2d& v);
把
struct Foo
{
Eigen::Vector2d v;
};
void my_function(Foo v);
改成
void my_function(const Foo& v);
编译器对堆栈对齐做出错误假设
局部解决方法:在定义函数的时候加上下面内容
__attribute__((force_align_arg_pointer)) void foo()
{
Eigen::Quaternionf q;
//...
}
全局解决方法:编译时添加-mincoming-stack-boundary=2
Pingback: C++库Eigen学习2:普通线性问题和分解 - Fivyex's Blog
Pingback: 机器人位姿反馈控制开发日志 - Fivyex's Blog