1
/* Ext. Eigen中的运算时机
2
  为提高效率,Eigen采用"懒求值"模式,也即使用各类运算符时,运算不会立即发生；
3
  相反,Eigen会返回一个表达式对象,实际计算在值实际被用到时才发生.
4
  "用到"值的最典型场景是将运算结果赋值给一个对象,也即使用operator=.
5
  有时会导致数据竞争,具体运算符条目中有说明.
6
*/
7
MatrixXi mat(3, 3);
8
mat << 1, 2, 3, 4, 5, 6, 7, 8, 9;
9
// 危险!对矩阵(1,1)点同时存在读取与写入操作,存在数据竞争.
10
mat.bottomRightCorner(2, 2) = mat.topLeftCorner(2, 2);
11
// 使用eval()函数显式生成中间变量,消除竞争.
12
mat.bottomRightCorner(2, 2) = mat.topLeftCorner(2, 2).eval();
13

14
// 1. 加法与减法
15
// 相关运算符: 单目运算符-; 双目运算符+ -; 复合运算符+= -=.
16
// 前提1:两对象的行数与列数对应
17
// 前提2:两对象的元素类型相同(自动类型提升不可用)
18
// 前提3:两侧运算对象均为array或均为matrix
19
Matrix2d m1;
20
Matrix3d m2;
21
Matrix2f m3;
22
Array22d a1, a2;
23
auto res = m1 + m2;  // 不合法,违反原则1
24
auto res = m1 + m3;  // 不合法,违反原则2
25
auto res = m1 + a1;  // 不合法,违反原则3
26
// 以下运算正确.
27
Matrix2d m4, m5;
28
m4 + m5;
29
m4 += m5;
30
m4 - m5;
31
m4 -= m5;
32
-m4;
33
m4 + a1.matrix();
34

35
// 2. 数乘与数除
36
// 以下运算均合法.
37
Matrix m1;
38
double d;
39
m1 * d;
40
d * m1;
41
m1 *= d;
42
m1 / d;
43
m1 /= d;
44

45
// 3. 转置与共轭
46
// 注意懒求值:a=a.transpose()会带来数据竞争,计算结果可能与预期不符.
47
Matrix3d m;
48
m.transpose();       // 转置(有数据竞争)
49
m.transposeInPlace();  // 原位转置(无数据竞争)
50
m.conjugate();       // 共轭
51
m.adjoint();         // 共轭转置(有数据竞争)
52
m.adjointInPlace();    // 原位共轭转置(无数据竞争)
53

54
// 4. 矩阵乘法
55
// 这是Eigen中唯一一种默认数据竞争会发生的情况.
56
// 如果实际上不会,可以在左侧添加.noalias()以提高性能.
57
Matrix3d m1, m2, m3;
58
m1 * m2;
59
m1 *= m2;  // == (m1 = m1 * m2)
60
m3.noalias() += m1 * m2;  // 如果确保不会发生数据竞争
61

62
// 5. 点乘与叉乘
63
Eigen::Vector3d v(1, 2, 3);
64
Eigen::Vector3d w(0, 1, 2);
65
v.dot(w);
66
v.cross(w);  // 对2维向量也可用,用以计算行列式
67

68
// 6. 基本算数操作
69
Eigen::Matrix3d m, m2;
70
Eigen::Vector3d v;
71
std::ptrdiff_t i, j;
72
m.sum();          // 所有元素的和
73
m.prod();         // 所有元素的积
74
m.mean();         // 所有元素的均值
75
m.minCoeff(&i, &j);    // 所有元素中的最小值(i,j是可选入参,表示位置)
76
v.minCoeff(&i);      // i是可选入参,表示位置
77
m.maxCoeff();        // 所有元素中的最大值
78
m.trace();        // 矩阵的迹
79
m.cwiseProduct(m2);   // m与m2的逐元素乘法,等效于(m.array() * m2.array()).matrix()
80
m.squaredNorm();    // 模的平方
81
m.norm();        // 模
82
m.normalize();      // m原位归一化
83
m.normalized();      // m非原位归一化
84
m.lpNorm<n>();      // n范数,(((m所有元素的绝对值)的n次方之和)的(1/n)次方.比如Norm是2范数.
85
m.lpNorm<Eigen::infinity>();  // +inf范数,m绝对值最大的元素.
86
m.sin();        // 逐元素sin(弧度制输入)
87
m.cos();
88
m.all();        // 返回true如果m元素全为真.
89
(a > 0).all();      // 合法.
90
m.any();        // 返回true如果m元素存在真.
91
m.count();        // 返回真元素数量.
92

93
// 7.向array转换,无运行时开销.
94
auto a = m.array();

1
// Array的大部分运算符与Matrix通用
2
// 区别在于Array间相乘时,返回逐元素积而非矩阵积.
3
Eigen::ArrayXXf a(2, 2);
4
Eigen::ArrayXXf b(2, 2);
5
a << 1, 2, 3, 4;
6
b << 5, 6, 7, 8;
7
auto c = a * b;  // c == {{5, 12}, {21, 32}};
8

9
// 2. Array的基本算数操作
10
Eigen::Array33d a;
11
a.abs();  // 返回|a|
12
a.sqrt(); // 逐元素开根号
13
a.min(a2);// 每个元素是a和a2中较小的一个
14

15
// 3.向Matrix转换,无运行时开销.
16
auto m = a.matrix();

1
// 视返回值的动/静态类型,分块方法略有不同.
2
// 注意:两种方法控制的是返回值类型,m矩阵本身动态/静态均可.
3
// 分块可以用作左值或右值.
4
Eigen::matrix4d m;
5
Eigen::vector4d v;
6
auto m2 = m.block(1,1,2,2); // m2是matrixXXd类型,以(1,1)为左上角取2x2的块
7
auto m3 = m.block<2,2>(1,1);// 分块大小编译期确定;m3是matrix22d类型,内容与m2一致.
8
m.block<2,2>(0,0) = m2;  // 合法.
9
m.block(0,0,2,2) = m3;   // 合法.
10
v.head(n);         // 矢量的头n个元素.
11
v.head<n>();
12
v.tail(n);         // 尾n个元素.
13
v.tail<n>();
14
v.segment(i, n);     // 从i开始的n个元素
15
v.segment<n>(i);
16

17
// 2.特殊分块函数
18
// 按行/按列/按角落分,比使用block方便高效.
19
// 2.1. 按行列分
20
// 返回值动/静态类型视m矩阵的动/静态类型.
21
Eigen::matrix4d m;
22
Eigen::matrixXd md;
23
m.row(i);   // 第i行,从0开始.静态矩阵.
24
m.col(j);   // 第j列,从0开始.静态矩阵.
25
md.row(i);  // 内容与上面相同.但是返回值类型为动态矩阵.
26
md.col(j);
27
m.topRows(q);     // 最上面q行
28
m.topRows<q>();
29
m.bottomRows(q);  // 最下面q行
30
m.bottomRows<q>();
31
m.leftCols(p);    // 左侧p列
32
m.leftCols<p>();
33
m.rightCols(p);    // 右侧p列
34
m.rightCols<p>();
35
m.middleCols(i, q); // 从i列开始的q列
36
m.middleCols<q>(i);
37
m.middleRows(i, p); // 从i行开始的q行
38
m.middleRows<p>(i);
39

40
// 2.2. 按角落分
41
Eigen::matrix4d m;
42
m.topLeftCorner(p, q);
43
m.topLeftCorner<p, q>();
44
m.bottomLeftCorner(p, q);
45
m.bottomLeftCorner<p, q>();
46
m.topRightCorner(p, q);
47
m.topRightCorner<p, q>();
48
m.bottomRightCorner(p, q);
49
m.bottomRightCorner<p, q>();
50

51
// 2.3. 对角
52
m.diagonal();  // (只适用方阵)取主对角元素,其余置为0.
53

54
// 3. Slicing, ()运算符的高级应用
55
/* 除了取单个元素, m(p,q)中的p,q还可以接受多种类型的参数.
56
 * 本质要求:定义了<integral type> operator[](<integral type>) const;
57
          <integral type> size() const; 两个函数.
58
 * 其中<integral type>可以是整型或者Eigen::Index可接受的类型(如std::ptrdiff_t).
59
 * 常见形式如下:
60
 * int:最常见的形式,取1行/1列.
61
 * Eigen::placeholders::all: 取该矩阵所有行/所有列.
62
 * ArithmeticSequence: 由Eigen::seq, Eigen::seqN或Eigen::lastN生成的序列.
63
 * Eigen数据结构:Eigen::Vector, Eigen::Matrix(仅限1维), Eigen::Array(仅限1维).
64
 * 标准库/C风格数据结构: std::Array, std::vector, int[N].
65
 */
66
// 3.1. seq与seqN的定义.
67
seq(firstIdx, lastIdx);      // seq(2, 6) = {2, 3, 4, 5, 6}
68
seq(firstIdx, lastIdx, incr);  // seq(2, 6, 2) = {2, 4, 6}
69
seqN(firstIdx, size);      // seqN(2, 4) = {2, 3, 4, 5}
70
seqN(firstIdx, size, incr);    // seqN(2, 4, 2) = {2, 4, 6, 8}
71
// 3.2. 取末位:lastN函数
72
A(lastN(p), lastN(q));  // == A.bottomRightCorner(p, q);
73
A(all, lastN(n,3));     // 取所有行,从最后一列开始,每隔2列取1列.
74
// 3.3. 编译期尺寸计算: 可以结合Eigen::fix()实现,在seq, seqN, lastN函数中.
75
v(seq(last-fix<7>, last-fix<2>);
76
v(seqN(last-7, fix<6>));
77
A(Eigen::placeholders::all, seq(fix<0>,last,fix<2>));  // A的所有偶数列(第一列为0列)
78
// 3.4. 倒序: 使用.reverse()或者seq/seqN.
79
A(seqN(last, n, fix<-1>), all);
80
A(lastN(n).reverse(), all);
81
// 3.5. 自定义下标列表
82
// 只需支持以下函数
83
/* Index s = ind.size(); or Index s = size(ind);
84
 * Index i;
85
 * i = ind[i];
86
 */
87
// 示例程序:pad
88
using Eigen::Index;
89
struct pad {
90
  Index size() const { return out_size; }
91
  Index operator[](Index i) const { return std::max<Index>(0, i - (out_size - in_size)); }
92
  Index in_size, out_size;
93
};
94

95
Matrix3i A;
96
A.reshaped() = VectorXi::LinSpaced(9, 1, 9);
97
// reshaped还没讲,此时A为{{1, 4, 7} ,{2, 5, 8}, {3, 6, 9}}.
98
MatrixXi B(5, 5);
99
B = A(pad{3, 5}, pad{3, 5});
100
/* B =
101
 * 1 1 1 4 7
102
 * 1 1 1 4 7
103
 * 1 1 1 4 7
104
 * 2 2 2 5 8
105
 * 3 3 3 6 9
106
*/

1
// 1. 预定义的高级Matrix/Array.
2
// 三种初始化方式.
3
Matrix3d::Zero();      // 0参数版本,适用于固定大小Matrix/Array.
4
ArrayXd::Zero(2);      // 1参数版本,适用于一维动态Matrix(Vector)/一维动态Array.
5
ArrayXXd::Zero(row, col);    // 2参数版本
6
// 若未专门说明,三种方式均可以用于以下方法.
7
// 方法列表.
8
ArrayXXd::Zero(row, col);    // 全0矩阵
9
ArrayXXd::Constant(row, col, val);  // 初值矩阵
10
ArrayXXd::Random(row, col);  // [0,1)随机初值
11
ArrayXXd::Identity(row, col);// 单位矩阵
12
VectorXd::Linspaced(size, low, high); // 线性排布,仅支持一维结构(size==1则返回high)
13
// 原位变换版本.
14
.setZero();
15
.setIdentity();
16
.setLinspaced(size, low, high);
17
// .finished()方法:基于"<<"运算符构建运算对象.
18
mat = (MatrixXf(2, 2) << 0, 1, 1, 0).finished() * mat;
19
// 2. 预定义的别名.
20
Vector3d::UnitX();
21
Vector3d::UnitY();
22
Vector3d::UnitZ();

1
// 1. 部分降阶.
2
// 1.1. 按行/列降阶.
3
.colwise()  // 按列降阶处理
4
.rowwise()  // 按行降阶处理
5
// 1.2. 常见用法
6
Eigen::Matrix3d m;
7
m.colwise().maxCoeff();  // 返回1*3向量,表示每列最大元素.
8
mat.colwise().sum().maxCoeff(&maxIndex); // 返回最大和列的序号
9

10
// 2. 广播.
11
// 第二操作数必须是Vector(对Matrix)或一维Array(对Array).
12
Eigen::MatriXd m(2, 4);
13
Eigen::VectorXd v(2);
14
mat.colwise() += v;  // 将v加到m的每一列.
15
(m.colwise() - v).colwise().squaredNorm().minCoeff(&index);  // 综合运用.

1
// Ex. 行优先与列优先
2

3
// Eigen 3.4: 添加reshaped(NRowsType,NColsType)函数,允许在原矩阵上创造视图.
4
// 提示:使用eval()以避免数据竞争.
5
Matrix4i m;
6
m.reshaped(2, 8);  // 按列取数,按列填充,(0,0)->(0,0),(1,0)->(1,0),(2,0)->(0,1), ...
7
m.reshaped();        // 默认形成单列向量
8
m.reshaped<RowMajor>();  // 按行取数,按列填充.
9
m.reshaped<AutoOrder>();// 根据默认行/列优先设置进行填充.
10

11
// 如果希望in-place改变存储的数据,使用resize().
12
Matrix4i m;
13
m.resize(2, 8);  // 根据默认行/列优先设置进行填充.

1
// 1D Vector/1D Array定义了随机访问迭代器.
2
// 节选支持的功能:
3
// Ranged-for
4
VectorXi v;
5
for (auto x : v) cout << x << " ";
6
// std::sort
7
std::sort(v.begin(), v.end());
8

9
// 多维数组可以使用reshaped()创建一维视图.
10
Matrix2i A;
11
for (auto x : A.reshaped()) cout << x << " ";
12

13
// 使用rowwise()/colwise()配合使用.
14
for (auto row : A.rowwise())
15
    std::sort(row.begin(), row.end());

1
/* 1. Map: 在C/C++数据结构上创建的Eigen视图
2
 * MapOptions: 内存对齐Aligned/Unaligned(默认).
3
 * StrideType<OuterStrideAtCompileTime_  ,
4
 *        InnerStrideAtCompileTime_>
5
 * 其中Inner表示相邻元素间的指针增量,Outer表示相邻行/列间的指针增量.
6
 * 两者均可指定为Dynamic以在运行时赋值.
7
 */
8
// 1.1. Map的基础定义
9
Map<Matrix<...> >
10
Map<Matrix<...>, int MapOptions, typename StrideType>
11
Map<const Matrix<...> >  // 只读Map
12
// 1.2. 基础用法
13
Map<MatrixXf> mf() = delete;  // 禁止默认初始化
14
Map<MatrixXf> mf(pf, rows, cols);  // pf:首地址float*类型指针
15
Map<const Vector4f> mi(pf);
16
// 1.3. Stride的使用示例.
17
// Stride<Dynamic, 2>(8, 2)表示Outer值为8,Inner值为2.
18
// 此处采用列优先存储,故同一列相邻行间的元素距离为Inner,
19
// 同一行相邻列之间的元素距离为Outer.
20
int array[24];
21
for (int i = 0; i < 24; ++i) array[i] = i;
22
cout << Map<MatrixXi, 0, Stride<Dynamic, 2> >(array, 3, 3, Stride<Dynamic, 2>(8, 2)) << endl;
23
/* res =
24
 *  0  8 16
25
 *  2 10 18
26
 *  4 12 20
27
 */
28
// 1.4. 修改Map对象指向的内存地址.
29
// 使用C++ Placement new语法.
30
int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
31
Map<RowVectorXi> v(NULL);
32
new (&v) Map<RowVectorXi>(data + 4, 5);  // 指向data+4位置长度为5的行向量.

1
// 1. 使用模板函数.
2
// 1.1. Eigen的所有运算对象都是少数基模板类的派生.
3
Eigen::MatrixBase  // 所有稠密矩阵运算对象的基类,用于仅矩阵函数.
4
Eigen::ArrayBase   // 所有稠密数组运算对象的基类,用于仅数组函数.
5
Eigen::DenseBase   // 上述两者的共同基类,用于适用于两者的函数.
6
Eigen::EigenBase   // 进一步包含特殊矩阵,稀疏矩阵等所有运算对象.
7

8
// 1.2. 示例函数
9
// 输出给定矩阵的逆条件数.
10
template <typename Derived>
11
void print_inv_cond(const MatrixBase<Derived>& a)
12
{
13
  const typename JacobiSVD<typename Derived::PlainObject>::SingularValuesType&
14
    sing_vals = a.jacobiSvd().singularValues();
15
  std::cout << "inv cond: " << sing_vals(sing_vals.size()-1) / sing_vals(0) << std::endl;
16
}
17
// 以下函数获取Array对象的最大元素.
18
template <typename Derived>
19
void print_max_coeff(const ArrayBase<Derived> &a)
20
{
21
  std::cout << "max: " << a.maxCoeff() << std::endl;
22
}
23
// 以下函数获取稠密运算对象的矩阵块.
24
template <typename Derived>
25
void print_block(const DenseBase<Derived>& b, int x, int y, int r, int c)
26
{
27
  std::cout << "block: " << b.block(x,y,r,c) << std::endl;
28
}
29
// 以下函数可以获取对象b的行列与总元素数量信息.
30
template <typename Derived>
31
void print_size(const Eigen::EigenBase<Derived>& b) {
32
  std::cout << "size (rows, cols): " << b.size() << " (" << b.rows() << ", " << b.cols() << ")" << std::endl;
33
}
34
// 多输入参数情况,需要给定不同的DeriveType.
35
template <typename DerivedA,typename DerivedB>
36
typename DerivedA::Scalar squaredist(const MatrixBase<DerivedA>& p1,const MatrixBase<DerivedB>& p2)
37
{
38
  return (p1-p2).squaredNorm();
39
}
40

41
// 2. 基于Eigen::ref非模板函数.
42
// 以下函数可以接受MatrixXf或其子块为实参.
43
float inv_cond(const Eigen::Ref<const Eigen::MatrixXf>& a) {
44
  const Eigen::VectorXf sing_vals = a.jacobiSvd().singularValues();
45
  return sing_vals(sing_vals.size() - 1) / sing_vals(0);
46
}
47

48
// 3. 在函数中resize数据对象
49
// 参考https://libeigen.gitlab.io/eigen/docs-nightly/TopicFunctionTakingEigenTypes.html,使用const_cast去掉const引用,过于黑魔法了,最好别用.