Eigen：评估从常量矩阵到新矩阵/数组的运算

Question

假设我有一个函数根据 Eigen 的文档模板化，以便使用

pybind11

在 C++ 和 Python 中使用它。

该函数的主要目标是执行笛卡尔 -> 极坐标转换。这是该函数的代码：

<typename VertexType>
typename Eigen::MatrixBase<VertexType>::PlainObject cartesian_to_polar(
    const Eigen::MatrixBase<VertexType>& cartesian
) {
    using MatrixType = typename Eigen::MatrixBase<VertexType>::PlainObject;
    using ScalarType = typename MatrixType::Scalar;
    using ReturnMatrixType = Eigen::Matrix<ScalarType, Eigen::Dynamic, Eigen::Dynamic, MatrixType::Options>;

    ReturnMatrixType polar;
    polar.resize(cartesian.rows(), 3);

    /* (1) */ const auto z = cartesian.col(2).eval().array();
    /* (2) */ const auto yx = (cartesian.col(1) / cartesian.col(0)).eval().array();

    polar.col(0) = Eigen::acos(z);
    polar.col(1) = Eigen::atan(yx);
    /* snip... compute radius and store in polar.col(2) */

    return polar;
}

根据 Eigen 关于惰性求值的规则，第 (1) 行从第三个

cartesian

矩阵的列中获取数据，然后

eval()

将其转换为一个新矩阵，可能是

MatrixBase<Scalar, Eigen::Dynamic, 1>

类型。但是，我在第 (2) 行遇到编译错误，其中错误消息指出我无法在两个

operator/

之间应用

ConstColXpr

。为什么？难道表达式的类型不应该抛弃它的常量性吗，因为表达式只能从源中读取？

我知道我无法修改原始（

cartesian

）矩阵，但是表达式的目标不是对矩阵应用符号运算吗？在这种情况下，原始块矩阵的常量性不应阻止我将一系列表达式计算为临时矩阵并对其进行修改。

---

而且，顺便说一句：我还在

Eigen::{acos,atan}(...)

行上遇到编译错误。如何将

MatrixBase

表达式显式转换为

ArrayBase

？当从编译时已知的

.array()

或依赖于 Python 的 Matrix<...>

 变量使用该函数时，

pybind11::EigenDRef<...> 方法似乎没有改变任何内容。

---

根据要求，这是使用该函数的最小示例：

#include "./cartesian_to_polar.h" /* Where the function above is defined. */
#include <Eigen/Eigen>
int main() {
    Eigen::Matrix3f vertices;
    vertices << 1.0f, .0f, .0f, .0f, 1.0f, .0f, .0f, .0f, 1.0f;
    auto polar_vertices = cartesian_to_polar(vertices);
}

Answer 1

首先，Eigen 和 C++11/14 的

auto

混合得不太好。这是因为 Eigen 在几乎所有情况下都会在调用方法时返回特殊的表达式模板，并且存储这些表达式模板可能会导致未定义的行为（例如，当临时变量的生命周期未延长时）。如果你知道自己在做什么，它就会起作用，但你应该小心。在您的具体情况下，这不是问题的根源，并且您在此处具体使用

auto

是安全的，但我强烈建议不要将

auto

与 Eigen 的 API 结合使用，但始终指定您想要的确切类型使用。 （您示例中的

auto

main()

没问题，因为您的 API 保证返回

PlainObject

。）

其次，您的问题与您模板化代码的事实无关。线

const auto yx = (cartesian.col(1) / cartesian.col(0)).eval().array();

失败，因为

cartesian.col(1)

和

cartesian.col(0)

是列向量表达式，而不是一维数组表达式，因此无法将它们相除。 （同样的方式，你也可以不做

VectorXd{...} / VectorXd{...}

。）你想要做的是在除法本身之前使用

.array()

：

const auto yx = (cartesian.array().col(1) / cartesian.array().col(0)).eval();

最后，在计算角度时，不应该使用

atan()

；相反，无论输入变量的符号如何，

atan2()

函数总是生成正确的值。 （例如，

atan(y/x)

会产生 y 和 x 负数的错误角度。）不幸的是，Eigen 没有为该函数提供预定义的包装器，因此您必须使用标准库：

polar.col(1) = cartesian.array().col(1).binaryExpr(cartesian.array().col(0), [] (ScalarType y, ScalarType x) { return std::atan2(y, x); });

而且不幸的是，

atan2()

通常比除法和普通

atan()

慢，但除非你绝对确定你的角度始终在预定的象限内，否则你仍然需要使用

atan2()

来获得正确的结果.

顺便说一句。你为什么要计算

acos(z)

？您是否尝试使用圆柱坐标或球坐标？无论哪种情况，

acos(z)

都是不正确的。 （在圆柱形中，您只需保持

z

不变，而在球形中，您可以计算

acos(z/r)

，其中

r

是 3D 中的总半径。）坐标的顺序也很奇怪 - 我个人会对它们进行排序规范地为

(rho, phi, z)

代替

(z, phi, rho)

（圆柱形情况），或

(r, theta, phi)

代替

(phi, theta, r)

（球形情况）。