Eigen3 教程基础篇（三）

参考

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;
  //...
}

这个就不测试了...

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