Eigen模板函数编译错误排查：坐标转换与延迟求值问题

我按照Eigen文档编写了一个兼容C++与pybind11的模板函数，用于实现笛卡尔坐标到极坐标的转换，代码如下：

template <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第三列求值为新矩阵，但第(2)行报错，提示无法在两个ConstColXpr间使用operator/。我有两个疑问：

1. 表达式仅读取源数据，为何const属性会阻碍运算？
2. Eigen::{acos,atan}行也存在编译错误，如何将MatrixBase表达式显式转为ArrayBase？.array()方法在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);
}

1. ConstColXpr无法使用operator/的原因 #

Eigen中，MatrixBase::col()返回的ConstColXpr（即Eigen::Block<const Derived>）是矩阵块的只读表达式，但报错和const属性无关——问题出在运算类型混淆：
Eigen的矩阵默认operator/是线性代数意义上的矩阵除法，两个列向量（矩阵块）之间做矩阵除法没有定义，哪怕是非const的ColXpr这么写也会报错。你需要的是逐元素除法，必须先将矩阵表达式转为Array表达式再运算。

正确写法应该先转Array再做除法：

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

若无需立即求值（Eigen延迟求值会自动优化），可省略.eval()：

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

2. 显式转换为ArrayBase并调用acos/atan的正确方式 #

Eigen的Eigen::acos/Eigen::atan是针对Array表达式的逐元素函数，不能直接传入Matrix表达式。你遇到的.array()无效问题，是调用时机错误：

• 对cartesian.col(2)这类矩阵块表达式，需先调用.array()转为Array表达式，再传入逐元素函数；
• pybind11的EigenDRef本质是Eigen矩阵的引用包装，完全兼容.array()/.matrix()接口，不存在失效问题。

修正后的完整代码：

template <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);

    // 先转Array视图，统一用逐元素运算
    const auto x_array = cartesian.col(0).array();
    const auto y_array = cartesian.col(1).array();
    const auto z_array = cartesian.col(2).array();

    // 对Array表达式调用逐元素函数，再转Matrix赋值给polar的列
    polar.col(0) = z_array.acos().matrix();
    polar.col(1) = (y_array / x_array).atan().matrix();
    // 计算半径：逐元素平方和开根号
    polar.col(2) = (x_array.square() + y_array.square() + z_array.square()).sqrt().matrix();

    return polar;
}

关键注意事项 #

• 严格区分矩阵运算（线性代数意义，如矩阵乘法*）和逐元素运算（需通过.array()切换，对应逐元素*//）；
• 优先使用Array表达式的成员函数（如array.acos()），比全局Eigen::acos(array)更清晰；
• 赋值给Matrix对象时，需用.matrix()将Array表达式转回Matrix表达式；
• pybind11的EigenDRef<T>完全兼容Eigen原生API，无需额外适配即可调用.array()。

内容的提问来源于stack exchange，提问作者Thibault de Villèle