eigen基础操作_eigen colwise

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>

1
2
3
4
5
6
7

基本运算

支持相同维数矩阵加减

支持乘除标量数

乘法（*）是矩阵乘法，不是对应元素相乘

基本方法

求转置、共轭矩阵、伴随矩阵

m.transpose();
m.conjugate();
m.adjoint();

注意：
不能是
m=m.transpose();//因为会在转置运算结束之前就开始把结果写进a
所以必须是：
b=a.transpose(); 
1
2
3
4
5
6
7
8
9

其他用法

v.dot(w); //点乘
v.cross(w);//叉乘
1
2

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()

1
2
3
4
5
6
7
8
9
10

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
1
2
3
4
5
6

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);
1
2
3
4
5
6
7
8

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";
1
2
3
4
5
6

另外，Eigen提供了MatrixXf::Index ，

MatrixXf::Index rownum,colnum;//可以存储行列号
1

Array（数组）使用

如果你需要做线性代数运算，如矩阵乘法，那么你应该使用矩阵;如果需要进行系数方面的操作，那么应该使用数组。

Array<float,3,3>a,b;
```
a*b;//此处是a与b各个元素相乘

1
2
3
4

//array与Matrix不能混用，需要转换类型
.array();
.matrix();
1
2
3

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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

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; //可以赋值   
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

void testblock() {
	Matrix<float, 2, 2> m;
	m = Matrix<float, 2, 2>::Constant(12);
	cout << m << endl;    //输出2*2的都是12的矩阵
}
1
2
3
4
5

//取某行或某列操作
m.col(2) += 3 * m.col(0);
1
2

对于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;
}
1
2
3
4
5
6
7
8
9

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
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

赋值方法

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;
1
2
3
4
5
6
7
8
9
10
11

2.特殊赋值

Array33f::Zero();
MatrixXf::Ones(3,2);
MatrixXf::Constant(5,2,1.58);
MatrixXf mat = MatrixXf::Random(2, 3);
1
2
3
4

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;
1
2
3
4
5
6
7

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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

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;
}
1
2
3
4
5
6
7
8
9

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 
}
1
2
3
4
5
6
7
8
9
10

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
1
2
3
4
5
6
7
8
9
10
11
12
13

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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

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";
}
1
2
3
4
5
6
7

目前还不知道如何将Map与std::vector直接转换

矩阵内部重叠问题

Eigen中用aliasing描述这一问题

1.在以下例子中重叠没有问题：

mat.bottomRightCorner(2,2) = mat.topLeftCorner(2,2);  //相当于将矩阵自身左上角四个元素赋值给右下角
1

2.但是当a = a.transpose();时就会因为重叠出现非编译错误

所以如何解决aliasing问题？

我们可以使用eval()函数：

a = a.transpose().eval();//原理大概是预先取一块内存空间，暂存a.transpose()结果，然后再复制给等式左边
1

当然对于转置，Eigen提供了：

a.transposeInPlace();//在原位置直接实现转置
1

其他常用方法的重叠解决方法

总结

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一样
1
2
3
4
5
6
7
8
9
10
11

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;
1
2
3
4
5
6
7

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
	}
}
1
2
3
4
5
6
7

如何选择：一般选择行存储，则按行遍历更快；Eigen默认是列，里面很多算法都是列优先，所以我们在选择时最好还是按照列优先

线性代数和分解

上表都是矩阵分解的方法，第一列是类名，第二列是类成员函数。这些类都有solve函数以及inverse函数

//以下形式便于理解类型
ColPivHouseholderQR<Matrix3f> dec(A);
Vector3f x = dec.solve(b);
1
2
3

HouseholderQR<MatrixXf> qr(50,50);
MatrixXf A = MatrixXf::Random(50,50);
qr.compute(A); // no dynamic memory allocation
qr.inverse();
1
2
3
4

//一般用以下方法求解$$Ax=b$$

Matrix2f x = A.ldlt().solve(b);
1

一般用以下方法求$A^{-1}$

A.colPivHouseholderQr().inverse();
1

检验显示：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;  //特征向量
}
1
2
3
4
5
6
7
8
9
10
11

注意：此时的特征值 特征向量不是从大到小排列

EigenSolver<MatrixXf> solver(covMat);
featureValue = solver.pseudoEigenvalueMatrix().diagonal();
featureVector = solver.pseudoEigenvectors();
1
2
3

稀疏矩阵

稀疏矩阵（Sparse Martix）需要include以下头文件

几何模块

#include <Eigen/Geometry> 
1

AlignedBox

这一模块用于包围盒操作

Eigen::AlignedBox< Scalar, AmbientDim >//<float,3>
1

简单实例：

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
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

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])); // 排序后的特征向量
	}
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24