Eigen3 教程基础篇（三）

最新推荐文章于 2024-09-29 21:48:15 发布

jucat

最新推荐文章于 2024-09-29 21:48:15 发布

阅读量666

点赞数 8

分类专栏： C++ 文章标签：机器学习人工智能机器人

本文链接：https://blog.csdn.net/jucat/article/details/142371724

版权

C++ 专栏收录该内容

4 篇文章 0 订阅

订阅专栏

参考

Eigen3 主页，Eigen3 官网教程

矩阵的本质，通过多种矩阵的应用去感受矩阵本质

3Blue1Brown 的线性代数，用可视化方法来表现线性代数的特性，强推

如何理解复数和虚数，有动画方便理解复数的意义

Eigen3 教程基础篇（一）

Eigen3 教程基础篇（二）

块操作

分块矩阵

一个分块矩阵是将矩阵分割出较小的矩形矩阵，这些较小的矩阵就称为区块。分块矩阵的分割原则是以水平线和垂直线进行划分。分块矩阵中，位在同一行（列）分区的每个子矩阵，都拥有相同的行数（列数）。

分块矩阵加法：如果两个分块矩阵的每个子矩阵行列相同，那么矩阵相加等于各个子矩阵相加。

分块矩阵乘法：分块矩阵 A 有 m*r 个子矩阵, 分块矩阵 B 有 r*n 个子矩阵，相乘后的矩阵 C 有 m*n 个子矩阵。

这里的 $A_{i,k}*B_{k,j}$ 同样要满足矩阵乘法规则。

矩阵取块操作

对矩阵 matrix 的 [ i, j ] 位置，取大小 ( p, q ) 的块操作：

动态大小的块操作：

matrix.block(i,j,p,q);

固定大小的块操作：

matrix.block<p,q>(i,j);

行和列是特殊的块，可以通过 row() 和 col()。

取出来的矩阵块可以视作独立的矩阵类进行操作。

void TestBlock()
{
  Eigen::Matrix4f m;
  m << 1.2, -2.3, 3.4, -4.5,
       5.6, -6.7, 7.8, -8.9,
       9.0, -0.1, 1.2, -2.3,
       3.4, -4.5, 5.6, -6.7;
  Eigen::Matrix2f md1122 = m.block(1,1,2,2);
  Eigen::Matrix2f mf1122 = m.block<2,2>(1,1);
  ROS_INFO_STREAM("动态取块，矩阵 m 的 (1,1) 取 2*2 块:" << std::endl << md1122);
  ROS_INFO_STREAM("固定取块，矩阵 m 的 (1,1) 取 2*2 块:" << std::endl << mf1122);

  Eigen::Matrix3f n;
  n << 1.2, -2.3, 3.4,
       -4.5, 5.6, -6.7,
       7.8, -8.9, 9.0;
  Eigen::Matrix2f nd0022 = m.block(0,0,2,2);
  Eigen::Matrix2f mn22 = md1122 + nd0022;
  ROS_INFO_STREAM("动态取块，矩阵 n 的 (0,0) 取 2*2 块:" << std::endl << nd0022);
  ROS_INFO_STREAM("矩阵 m(1,1){2*2} + n(0,0){2*2}" << std::endl << mn22);

  Eigen::Matrix2f md2022 = m.block(2,0,2,2);
  Eigen::Matrix2f nf1122 = n.block<2,2>(1,1);
  ROS_INFO_STREAM("矩阵 m(2,0){2*2} * n(1,1){2*2}" << std::endl << md2022*nf1122);
}

向量取块操作

取 n 个向量头部元素的块，动态块操作：vector.head(n)，固定区块：vector.head<n>()；

取 n 个向量尾部元素的块，动态块操作：vector.tail(n)，固定区块：vector.tail<n>()；

在 i 位置取 n 个元素的块，动态块操作：vector.segment(i,n)，固定区块：vector.segment<n>(i)；

void TestVectorBlock()
{
  Eigen::Vector3f v3f(1.2, -2.3, 3.4);
  Eigen::Vector2f v3fh2 = v3f.head(2);
  Eigen::Vector2f v3ft2 = v3f.tail(2);
  ROS_INFO_STREAM("向量 v3f:" << std::endl << v3f);
  ROS_INFO_STREAM("向量 v3f 取头部 2 个元素:" << std::endl << v3fh2);
  ROS_INFO_STREAM("向量 v3f 取尾部 2 个元素:" << std::endl << v3ft2);

  Eigen::VectorXf v6f(6);
  v6f << 0.1, -1.2, 2.3, -3.4, 4.5, -5.6;
  Eigen::Vector3f v6f33 = v6f.segment(3,3);
  Eigen::Vector3f v6f23 = v6f.segment<3>(2);
  ROS_INFO_STREAM("向量 v6f:" << std::endl << v6f);
  ROS_INFO_STREAM("向量 v6f 动态取块在 3 位置取 3 个元素:" << std::endl << v6f33);
  ROS_INFO_STREAM("向量 v6f 在 2 位置固定取块 3 个元素:" << std::endl << v6f23);
}

Reductions 归约

在 Eigen 中，Reductions（归约）是返回单个标量值的函数方法。

常见的归约方法：

sum() 矩阵系数和，prod() 矩阵系数乘积，mean() 矩阵系数均值，minCoeff() 矩阵中最小系数，maxCoeff() 矩阵中最大系数，trace() 矩阵中对角线（左上角到右下角）元素之和。

squaredNorm() 向量系数绝对值的平方之和，norm() 向量系数绝对值平方和的平方根。

这两个方法对矩阵同样适用，此时矩阵 M(n,p) 被看做 n*p 大小的向量。

布尔归约：

all()，如果矩阵，数组的所有系数均计算为 true，返回 true；

any()，如果矩阵，数组的存在任一系数计算为 true，返回 true；

count()，返回矩阵或数组中系数计算为 true 的数量；

部分归约：

colwise() 和 rowwise() 可以取得列向，或行向的归约。

void TestReductions()
{
  Eigen::Matrix4f m;
  m << 1.2, -2.3, 3.4, -4.5,
       5.6, -6.7, 7.8, -8.9,
       9.0, -0.1, 1.2, -2.3,
       3.4, -4.5, 5.6, -6.7;
  ROS_INFO_STREAM("矩阵 m:" << std::endl << m);

  ROS_INFO_STREAM("m.sum():" << m.sum() << std::endl <<
                  "m.prod():" << m.prod() << std::endl <<
                  "m.mean():" << m.mean() << std::endl <<
                  "m.minCoeff():" << m.minCoeff() << std::endl <<
                  "m.maxCoeff():" << m.maxCoeff() << std::endl <<
                  "m.trace():" << m.trace() << std::endl <<
                  "m.squaredNorm():" << m.squaredNorm() << std::endl <<
                  "m.norm():" << m.norm() << std::endl);

  ROS_INFO_STREAM("m.colwise().maxCoeff():" << m.colwise().maxCoeff() << std::endl <<
                  "m.rowwise().sum().minCoeff():" << m.rowwise().sum().minCoeff() << std::endl);
}

broadcasting 传播

对所有的方向进行操作。

void TestBroadcasting()
{
  Eigen::MatrixXf m(2,4);
  m << 1.0, -11.23, 6.9, -0.83,
       -5.6, 2.56, -2.97, 6.56;
  ROS_INFO_STREAM("矩阵 m:" << std::endl << m);

  Eigen::VectorXf n(2);
  n << 2,
       3;
  ROS_INFO_STREAM("向量 n:" << std::endl << n);

  Eigen::Index index;
  (m.colwise() - n).colwise().squaredNorm().minCoeff(&index);
  std::cout << "Nearest neighbour is column " << index << ":" << std::endl;
  std::cout << m.col(index) << std::endl;
}

Map 类：原始缓存数据接口

当在 Eigen 之外定义了一组数据，这组数据就是一个矩阵或者向量。这里可以通过 Eigen 的 Map 模板来使用 Eigen 的算法处理这组数据。

Map 模板类：

Map<Matrix<typename Scalar, int RowsAtCompileTime, int ColsAtCompileTime> >

例如Map<MatrixXf> mf(pf,rows,cols)，pf 是 Eigen 之外定义的 rows*cols 大小的矩阵的数据内存地址，这里通过 Map 模板构造的对象 mf 得到了和 MatrixXf 等效的数据处理接口，可以像使用 MatrixXf 一样处理 pf 地址的数据。

修改 Map 模板类对象的数据地址，挺奇怪的写法：

new (&map) Map<Matrix<typename Scalar, int Rows, int Cols> >(new_data);

void TestMap()
{
  int array[8];
  for(int i = 0; i < 8; ++i)
    array[i] = i;
  auto m = Eigen::Map<Eigen::Matrix<int,2,4>>(array);
  ROS_INFO_STREAM("矩阵 array -> m:" << std::endl << m);
  ROS_INFO_STREAM("矩阵 array -> m.colwise().sum()" << std::endl << m.colwise().sum());

  int data[] = {1,2,3,4,5,6,7,8,9};
  new (&m) Eigen::Map<Eigen::Matrix<int,2,4>>(data);
  ROS_INFO_STREAM("矩阵 data -> m:" << std::endl << m);
  ROS_INFO_STREAM("矩阵 data -> m.rowwise().squaredNorm():" << std::endl << m.rowwise().squaredNorm());
}

Aliasing 混叠

在 Eigen 中，Aliasing 混叠指赋值语句左右两侧出现相同的矩阵/数组/向量，这里指左右操作数内存上有重叠。有时候 Aliasing 混叠会导致异常。

混叠异常情况

void TestAliasing()
{
  Eigen::MatrixXi m(3,3); 
  m << 1, 2, 3,
       4, 5, 6,
       7, 8, 9;
  ROS_INFO_STREAM("矩阵 m:" << std::endl << m);
  
  m.bottomRightCorner(2,2) = m.topLeftCorner(2,2);
  ROS_INFO_STREAM("after the assignment, m:" << std::endl << m);
}

上面的例子中，系数 5 同时存在于左上 2*2 矩阵块和右下 2*2 矩阵块中，Eigen 一般情况下为了性能直接修改数据内存，而没有拷贝右值的临时变量处理。这就导致了混叠操作有可能出现计算错误。

其他像 a = a.transpose()的操作也有导致出错！

避免混叠异常

通过 eval() 函数给右值创建临时对象。

void TestAliasing()
{
  ROS_INFO_STREAM("---------- 混叠异常 ----------");
  Eigen::MatrixXi m(3,3); 
  m << 1, 2, 3,
       4, 5, 6,
       7, 8, 9;
  auto x = m;
  ROS_INFO_STREAM("矩阵 m:" << std::endl << m);
  
  m.bottomRightCorner(2,2) = m.topLeftCorner(2,2);
  ROS_INFO_STREAM("after the assignment, m:" << std::endl << m);

  ROS_INFO_STREAM("---------- 避免混叠 ----------");
  x.bottomRightCorner(2,2) = x.topLeftCorner(2,2).eval();
  ROS_INFO_STREAM("after the assignment, x:" << std::endl << x);
}

另外还有 ***InPlace() 接口避免混叠异常：

矩阵乘法是唯一被默认假定混叠的运算，所以矩阵乘法是默认通过临时对象计算，此时混叠不会出现异常。

也是由于矩阵乘法的特殊处理，其中的临时对象赋值和销毁消耗了时间，此时可以通过 noalias() 函数指定矩阵乘法不使用临时对象。

下面用 noalias() 告诉 eigen 库，左值的矩阵数据内存不与右值操作数有重叠，不需要拷贝右值的临时变量。

void TestAliasingMultiplication()
{
  Eigen::MatrixXi m(3,3); 
  m << 1, 2, 3,
       4, 5, 6,
       7, 8, 9;
  auto x = m;
  ROS_INFO_STREAM("矩阵 m:" << std::endl << m);

  auto mx = m * x;
  ROS_INFO_STREAM("矩阵 m * 矩阵 x:" << std::endl << mx);

  m *= m;
  ROS_INFO_STREAM("矩阵 m 乘法混叠无异常, m:" << std::endl << m);

  x.noalias() = x * x;
  ROS_INFO_STREAM("矩阵 x 混叠，不使用临时变量, x:" << std::endl << x);
}

一般来说，可以用 noalias() 指定不创建临时对象，加速计算；用 eval() 来创建临时对象，避免计算错误。

Explanation of the assertion on unaligned arrays 未对齐内存问题

使用 eigen 库时候，可能出现下面的断言错误：

这个错误的根本原因是固定矩阵/向量大小的 Eigen 对象没有在正确对齐的位置创建，导致 SIMD 指令寻址错误。

！！！使用 C++17 可以避免这个问题！！！

！！！动态大小的矩阵不会有这个问题，Eigen 会管理动态内存的对象！！！

常见的出错场景（C++11编译测试）

构建具有固定大小的 Eigen 成员变量的对象

class xclass
{
 public:
  Eigen::Matrix2d v;
};

void TestUnalignedArray()
{
  auto x = new xclass();
  x->v << 1.0, 2.0, 3.0, 4.0;
  ROS_INFO_STREAM("x.v:" << std::endl << x->v);
}

未出现上面的错误：

使用容器存储 Eigen 对象

void TestUnalignedArray()
{
  std::vector<Eigen::Vector2d> ve2d;
  Eigen::Vector2d m2d0;
  m2d0 << 1.0, 2.0;
  Eigen::Vector2d m2d1;
  m2d1 << 5.0, 6.0;

  ve2d.push_back(m2d0);
  ve2d.push_back(m2d1);
  ROS_INFO_STREAM("ve2d.front()" << std::endl << ve2d.front());
  ROS_INFO_STREAM("ve2d.at(1)" << std::endl << ve2d.at(1));
}

未出现上面的错误：

按值传递 Eigen 对象

void PrintInputMatrix(const Eigen::Matrix2d m)
{
  ROS_INFO_STREAM("input m:" << std::endl << m);
}

void TestUnalignedArray()
{
  Eigen::Matrix2d m;
  m << 1.2, 2.3, 3.4, 4.5;
  PrintInputMatrix(m);
}

未出现上面的错误：

在堆栈创建 Eigen 对象（Win 中的 gcc）

void foo()
{
  Eigen::Quaternionf q;
  //...
}

这个就不测试了...

不知道为啥就是没出现异常 :) 具体的解决办法参考原文资料吧。

jucat

关注

8
点赞
踩
25

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录

Eigen3 教程基础篇（三）

参考

相关文章

块操作

分块矩阵

矩阵取块操作

向量取块操作

Reductions 归约

broadcasting 传播

Map 类：原始缓存数据接口

Aliasing 混叠

Explanation of the assertion on unaligned arrays 未对齐内存问题

常见的出错场景（C++11编译测试）