在 Eigen C++ 中广播行和列向量

Question

我在 NumPy 中编写了以下 Python 代码：

> r = 3
> y, x = numpy.ogrid[-r : r + 1, -r : r + 1]
> mask = numpy.sqrt(x**2 + y**2)
> mask

array([[4.24264, 3.60555, 3.16228, 3.00000, 3.16228, 3.60555, 4.24264],
       [3.60555, 2.82843, 2.23607, 2.00000, 2.23607, 2.82843, 3.60555],
       [3.16228, 2.23607, 1.41421, 1.00000, 1.41421, 2.23607, 3.16228],
       [3.00000, 2.00000, 1.00000, 0.00000, 1.00000, 2.00000, 3.00000],
       [3.16228, 2.23607, 1.41421, 1.00000, 1.41421, 2.23607, 3.16228],
       [3.60555, 2.82843, 2.23607, 2.00000, 2.23607, 2.82843, 3.60555],
       [4.24264, 3.60555, 3.16228, 3.00000, 3.16228, 3.60555, 4.24264]])

现在，我正在 Eigen 中制作掩码，我需要在其中广播行向量和列向量。不幸的是，这是不允许的，所以我做了以下解决方法：

int len = 1 + 2 * r;
MatrixXf mask = MatrixXf::Zero(len, len);
ArrayXf squared_yx = ArrayXf::LinSpaced(len, -r, r).square();
mask = (mask.array().colwise() + squared_yx) + 
       (mask.array().rowwise() + squared_yx.transpose());
mask = mask.cwiseSqrt();
cout << "mask" << endl << mask << endl;

4.24264 3.60555 3.16228       3 3.16228 3.60555 4.24264
3.60555 2.82843 2.23607       2 2.23607 2.82843 3.60555
3.16228 2.23607 1.41421       1 1.41421 2.23607 3.16228
      3       2       1       0       1       2       3
3.16228 2.23607 1.41421       1 1.41421 2.23607 3.16228
3.60555 2.82843 2.23607       2 2.23607 2.82843 3.60555
4.24264 3.60555 3.16228       3 3.16228 3.60555 4.24264

有效。但我想知道是否有另一种更短的方法来做到这一点。因此我的问题是如何在 Eigen C++ 中广播行和列向量?

系统信息

Tool	Version
Eigen	3.3.7
GCC	9.4.0
Ubuntu	20.04.4 LTS

Answer 1

我认为最简单的方法（如：最具可读性）是复制。

  int r = 3;
  int len = 1 + 2 * r;
  const auto& squared_yx = Eigen::ArrayXf::LinSpaced(len, -r, r).square();
  const auto& bcast = squared_yx.replicate(1, len);
  Eigen::MatrixXf mask = (bcast + bcast.transpose()).sqrt();

请注意，您所做的事情在数值上是不稳定的（对于较大的 r）并且存在 hypot 函数来解决这些问题。所以即使你的 python 代码也可以更好：

r = 3
y, x = numpy.ogrid[-r : r + 1, -r : r + 1]
mask = numpy.hypot(x, y)

要在 Eigen 中实现相同的效果，请执行以下操作：

  const auto& yx = Eigen::ArrayXf::LinSpaced(len, -r, r);
  const auto& bcast = yx.replicate(1, len);
  Eigen::MatrixXf mask = bcast.binaryExpr(bcast.transpose(),
      [](float x, float y) noexcept -> float {
          return std::hypot(x, y);
  });

Eigen 关于 binaryExpr 的文档目前已损坏，因此很难找到。

公平地说，在这种特殊情况下，您可能永远不会运行陷入稳定性问题，因为您首先会运行内存不足。但是，它仍然想指出这一点，因为看到天真的 sqrt(x**2 + y**2) 总是有点危险。此外，在 Python 中，hypot 从性能角度来看可能仍然值得，因为它减少了临时内存分配和函数调用的数量。

二进制表达式

缺少关于 binaryExpr 的文档，我猜是因为解析器在处理 Eigen 的 C++ 代码时遇到了问题。在任何情况下，都可以间接找到它 CwiseBinaryOp and similarly CwiseUnaryOp, CwiseNullaryOp and CwiseTernaryOp.

使用看起来有点奇怪但是很简单。它需要一个仿函数（一个带有 operator() 的结构、一个函数指针或一个 lambda）并应用这个 element-wise.

一元运算使这一点非常清楚。如果 Eigen::Array.sin() 不存在，你可以这样写： array.unaryExpr([](double x) -> double { return std::sin(x); })达到完全一样的效果。

二元和三元版本将一个或两个以上的本征表达式作为函数的第二个和第三个参数。这就是我上面所做的。 nullary 版本在文档 in its own chapter.

中进行了说明

使用自动

Eigen 警告 auto 是正确的，但前提是您必须知道自己在做什么。重要的是要意识到 Eigen 表达式上的 auto 只是使表达式保持不变。它不会将其计算为向量或矩阵。

如果您想编写一个在单个语句中难以阅读的复杂表达式，这很好而且非常有用。在我上面的代码中，没有临时内存分配，也没有进行浮点计算，直到将最终表达式分配给矩阵。

只要程序员知道这些是表达式而不是最终矩阵，一切都很好。

我认为主要的 take-away 是将 auto 与 Eigen 一起使用应限于 short-lived（如：在单个函数内）标量表达式。任何对一切都使用 auto 的编码风格将很快被 Eigen 破坏或难以阅读。但它可以安全地使用，并使代码在过程中更具可读性，而不会像评估矩阵那样牺牲性能。

至于为什么我选择const auto&而不是auto或const auto：主要是与手头任务无关的习惯力量。我主要针对这样的情况这样做：

const Eigen::Vector& Foo::get_bar();

void quz(Foo& foo)
{
    const auto& bar = foo.get_bar();
}

在这里，bar 将保留一个引用，而 auto 将创建一个副本。如果更改了 return 值，一切都将保持有效。

Eigen::Vector Foo::get_bar();

void quz(Foo& foo)
{
    const auto& bar = foo.get_bar();
}

现在无论如何都创建了一个副本。但是一切都会继续工作，因为将 return 值分配给 const-reference 会延长对象的生命周期。所以这可能看起来像一个悬挂指针，但它不是。

在 Eigen C++ 中广播行和列向量

Broadcasting Row and Column Vector in Eigen C++

c++

eigen

eigen3

array-broadcasting

系统信息

二进制表达式

使用自动