Matrix基本使用
Eigen概述 - Go吧LpengSu | Blog (yueyuebird-su.github.io)
Eigen 学习笔记(一)_逐梦者-CSDN博客_eigen
Eigen中Matrix 与Vector相似
Matrix<typename Scalar, //类型
int RowsAtCompileTime,
int ColsAtCompileTime,
int Options = 0, //默认是列优先
int MaxRowsAtCompileTime = RowsAtCompileTime,
int MaxColsAtCompileTime = ColsAtCompileTime>
基本运算
支持相同维数矩阵加减
支持乘除标量数
乘法(*)是矩阵乘法,不是对应元素相乘
基本方法
求转置、共轭矩阵、伴随矩阵
m.transpose();
m.conjugate();
m.adjoint();
注意:
不能是
m=m.transpose();//因为会在转置运算结束之前就开始把结果写进a
所以必须是:
b=a.transpose();
其他用法
v.dot(w); //点乘
v.cross(w);//叉乘
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;//矩阵的迹等价于mat对角化后对角线元素之和a.diagonal().sum()
Here is mat.sum(): 10 Here is mat.prod(): 24 Here is mat.mean(): 2.5 Here is mat.minCoeff(): 1 Here is mat.maxCoeff(): 4 Here is mat.trace(): 5
Eigen: Eigen::MatrixBase< Derived > Class Template Reference
访问元素
MatrixXf m(2, 2);
m = MatrixXf::Random(2, 2);
cout << m(1, 1) << endl;
cout << m[0][0] << endl; //这样是错误的,对于matrix必须是括号访问
VectorXd v(3);
v = Vector3d::Random(3);
cout << v[0];
cout << v(1);
Matrix3f m = Matrix3f::Random();
std::ptrdiff_t i, j;
float minOfM = m.minCoeff(&i,&j); //将m矩阵中最小数的位置传到i j ;把最小数传到minOfm
cout << "Here is the matrix m:\n" << m << endl;
cout << "Its minimum coefficient (" << minOfM
<< ") is at position (" << i << "," << j << ")\n\n";
另外,Eigen提供了MatrixXf::Index ,
MatrixXf::Index rownum,colnum;//可以存储行列号
Array(数组)使用
如果你需要做线性代数运算,如矩阵乘法,那么你应该使用矩阵;如果需要进行系数方面的操作,那么应该使用数组。
Array<float,3,3>a,b;
```
a*b;//此处是a与b各个元素相乘
//array与Matrix不能混用,需要转换类型
.array();
.matrix();
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;
//Eigen allows assigning array expressions to matrix variables
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;
}
block分块
int main()
{
Eigen::MatrixXf m(4,4),a(2,2);
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;
}
m.block<2,2>(1,1)=a; //可以赋值
}
Block in the middle 6 7 10 11 Block of size 1x1 1 Block of size 2x2 1 2 5 6 Block of size 3x3 1 2 3 5 6 7 9 10 11
void testblock() {
Matrix<float, 2, 2> m;
m = Matrix<float, 2, 2>::Constant(12);
cout << m << endl; //输出2*2的都是12的矩阵
}
//取某行或某列操作
m.col(2) += 3 * m.col(0);
对于Vector
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;
}
v.head(3) = 1 2 3 v.tail<3>() = 4 5 6 after 'v.segment(1,4) *= 2', v = 1 4 6 8 10 6
赋值方法
1.“<<”+逗号赋值法 前提:待赋值矩阵必须大小确定
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;
2.特殊赋值
Array33f::Zero();
MatrixXf::Ones(3,2);
MatrixXf::Constant(5,2,1.58);
MatrixXf mat = MatrixXf::Random(2, 3);
ArrayXXf table(10, 4);
table.col(0) = ArrayXf::LinSpaced(10, 0, 90); //0--90 取10个数
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;
Reductions, visitors and broadcasting
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; //计算平方范数 5
cout << "v.norm() = " << v.norm() << endl; //计算平方范数的根 2.23607
cout << "v.lpNorm<1>() = " << v.lpNorm<1>() << endl; //求取矩阵一范数 3
cout << "v.lpNorm<Infinity>() = " << v.lpNorm<Infinity>() << endl; //求取无穷范数 2
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;
}
void test03() {
ArrayXXf a(2, 2);
a << 1, 2,
3, 4;
cout << "(a > 0).all() = " << (a > 0).all() << endl; //1
cout << "(a > 0).any() = " << (a > 0).any() << endl; //1
cout << "(a > 0).count() = " << (a > 0).count() << endl; //4
cout << endl;
}
colwise() and broadcasting
void test04() {
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.maxCoeff();返回的是整个mat最大值,此时返回值只有9;
//但是有了colwise(),则为按列求取最大值,返回值为各列最大值 3 2 7 9
}
colwise还可以实现如下功能
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;
//1 2 6 9
//4 2 8 3
col-wise同理
需要注意的是:“+”后面的那个必须是Vector,否则程序报错
broadcasting与其他操作可以结合
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; //0
cout << m.col(index) << endl; //1 3
}
Map
Map中提供了Eigen与C++基本数组转换的方法,即通过指针与所占内存空间进行转换
void testMap02() {
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与std::vector直接转换
矩阵内部重叠问题
Eigen中用aliasing描述这一问题
1.在以下例子中重叠没有问题:
mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2); //相当于将矩阵自身左上角四个元素赋值给右下角
2.但是当a = a.transpose();
时就会因为重叠出现非编译错误
所以如何解决aliasing问题?
我们可以使用eval()函数:
a = a.transpose().eval();//原理大概是预先取一块内存空间,暂存a.transpose()结果,然后再复制给等式左边
当然对于转置,Eigen提供了:
a.transposeInPlace();//在原位置直接实现转置
其他常用方法的重叠解决方法
总结
1.alias对于标量计算、array和Matrrix加法是没有威胁的
2.当你使用矩阵乘法时,他默认会发生重叠。当你肯定重叠不会产生影响时,则可以使用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;
//两个B一样
3.其他情况下,Eigen都默认不会重叠,所以需要判断,视情况加.eval();
4.从版本3.3之后,对于矩阵乘法,如果矩阵shape改变并且乘积结果不是直接赋值给左边,则也默认不会重叠!!
MatrixXf A(2,2), B(3,2);
B << 2, 0, 0, 3, 1, 1;
A << 2, 0, 0, -2;
A = (B * A).cwiseAbs(); //由于B*A后面跟了方法,不是直接赋值,所以此处发生重叠,结果错误
//必须是:
A = (B * A).eval().cwiseAbs();
cout << A;
ColMajor&RowMajor
void ColRowMajor() {
Matrix<int, 3, 3, ColMajor> m;
m << 1, 2, 3, 4, 5, 6, 7, 8, 9;
for (int i = 0; i < 9; i++) {
cout << *(m.data() + i) << endl;//结果是:1 4 7 2 5 8 3 6 9
}
}
如何选择:一般选择行存储,则按行遍历更快;Eigen默认是列,里面很多算法都是列优先,所以我们在选择时最好还是按照列优先
线性代数和分解
上表都是矩阵分解的方法,第一列是类名,第二列是类成员函数。这些类都有solve函数以及inverse函数
//以下形式便于理解类型
ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);
HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation
qr.inverse();
//一般用以下方法求解 A x = b Ax=b Ax=b
Matrix2f x = A.ldlt().solve(b);
一般用以下方法求 A − 1 A^{-1} A−1
A.colPivHouseholderQr().inverse();
检验显示:m.colPivHouseholderQr().inverse()
求解逆矩阵比m.inverse()
更准确
特征值 特征向量
void eigenvalues() {
Matrix2f A;
A << 1, 2, 2, 3;
cout << "Here is the matrix A:\n" << A << endl;
SelfAdjointEigenSolver<Matrix2f> eigensolver(A);
if (eigensolver.info() != Success) abort();
cout << "The eigenvalues of A are:\n" << eigensolver.eigenvalues() << endl; //特征值
cout << "Here's a matrix whose columns are eigenvectors of A \n"
<< "corresponding to these eigenvalues:\n"
<< eigensolver.eigenvectors() << endl; //特征向量
}
注意:此时的特征值 特征向量不是从大到小排列
EigenSolver<MatrixXf> solver(covMat);
featureValue = solver.pseudoEigenvalueMatrix().diagonal();
featureVector = solver.pseudoEigenvectors();
稀疏矩阵
稀疏矩阵(Sparse Martix)需要include以下头文件
几何模块
#include <Eigen/Geometry>
AlignedBox
这一模块用于包围盒操作
Eigen::AlignedBox< Scalar, AmbientDim >//<float,3>
简单实例:
void testAlignedBox() {
AlignedBox2d myAlignedBox;
Vector2d m;
m << 1, 0.5;
myAlignedBox.extend(m);
m << 2, 2;
myAlignedBox.extend(m);
cout <<"max():"<<myAlignedBox.max() << endl;//返回最大的corner
m << 3, 4;
cout << "exteriorDistance():" << myAlignedBox.exteriorDistance(m) //外部一个向量m与包围盒距离
<< "squaredExteriorDistance():" << myAlignedBox.squaredExteriorDistance(m) << endl;//外部一向量与包围盒距离平方
//再构建一个alignedBox
AlignedBox2d myAlignedBox2,myalignedBox3;
Vector2d b;
b << 0.5, 0.5;
myAlignedBox2.extend(b);
b << 1.5, 1.5;
myAlignedBox2.extend(b);
myalignedBox3 = myAlignedBox2.intersection(myAlignedBox);//1 2 做交集 此外,intersects()返回的是bool值
cout << myalignedBox3.BottomLeft << endl;
b << 1.25, 1;
cout <<"contains (0 or 1):"<< myalignedBox3.contains(b) << endl;//看b是否在myalignedBox3中,参数也可以是box
cout << "diagonal():" << myalignedBox3.diagonal() << endl;//对角线向量
cout << "random sample:" << myAlignedBox.sample() << endl;//以均匀分布抽样的边界框内的随机点
cout << "sizes():" << myalignedBox3.sizes() << endl;//myalignBox3的长宽//vector
}
cornerType:
transform和translate稍后再看
AngleAxis
注意:AngleAxis运用时:Axis的轴需要归一化
首先,Eigen在MatrixBase中定义了
UnitX(); UnitY(); UnitZ(); UnitW(); 作为默认的四个轴
例如:
Vector3d::UnitX();
必须是固定大小的
附录
将特征值特征向量按特征值大小排列
using eigMatrix = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic>;
using eigVector = Eigen::Matrix<float, Eigen::Dynamic, 1>;
using stdTupleEigen = std::tuple<float, eigMatrix>;
using stdVectorTuple = std::vector<stdTupleEigen>;
void sortEigenVectorByValues(eigVector& eigenValues, eigMatrix& eigenVectors)
{
stdVectorTuple eigenValueAndVector;
int size = static_cast<int>(eigenValues.size()); //强制转换
eigenValueAndVector.reserve(size);
for (int i = 0; i < size; ++i)
eigenValueAndVector.push_back(stdTupleEigen(eigenValues[i], eigenVectors.col(i)));
// 使用标准库中的sort,按从大到小排序
std::sort(eigenValueAndVector.begin(), eigenValueAndVector.end(),
[&](const stdTupleEigen& a, const stdTupleEigen& b) -> bool {
return std::get<0>(a) > std::get<0>(b);
});
for (int i = 0; i < size; ++i) {
eigenValues[i] = std::get<0>(eigenValueAndVector[i]); // 排序后的特征值
eigenVectors.col(i).swap(std::get<1>(eigenValueAndVector[i])); // 排序后的特征向量
}
}