在Eigen中编写惰性求值无拷贝高性能函数的技术咨询

一、核心思路：用Eigen模板特性替代预处理器冗余 #

你之前靠预处理器为每个变量定义函数的方式，确实会让代码变得臃肿难维护。Eigen的模板元编程和表达式模板是天生的高性能工具，能在编译期完成优化，完全不需要重复定义那套操作。

二、类numpy.diff()的有限差分函数实现 #

先明确numpy.diff()的核心逻辑：沿指定轴计算相邻元素的差值（默认轴为最后一维）。我们用模板函数实现，同时适配不同矩阵类型和差分轴。

1. 基础1D行差分（对应numpy.diff(arr, axis=0)） #

#include <Eigen/Core>

template <typename Derived>
Eigen::Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime - 1, Derived::ColsAtCompileTime>
diff_1d(const Eigen::MatrixBase<Derived>& mat) {
    // 取第1行到末尾的元素，减去第0行到倒数第二行的元素
    return mat.bottomRows(mat.rows() - 1) - mat.topRows(mat.rows() - 1);
}

2. 支持指定轴的通用差分函数 #

如果要像numpy一样支持axis参数（0=行方向，1=列方向），用编译期分支实现零开销：

template <int Axis = 1, typename Derived>
auto diff(const Eigen::MatrixBase<Derived>& mat) {
    static_assert(Axis == 0 || Axis == 1, "Axis must be 0 (rows) or 1 (cols)");
    
    if constexpr (Axis == 0) {
        return mat.bottomRows(mat.rows() - 1) - mat.topRows(mat.rows() - 1);
    } else { // Axis == 1
        return mat.rightCols(mat.cols() - 1) - mat.leftCols(mat.cols() - 1);
    }
}

使用示例：

int main() {
    Eigen::MatrixXd mat(3, 3);
    mat << 1, 2, 3,
           4, 5, 6,
           7, 8, 9;
    
    // 行方向差分（axis=0）
    auto diff_rows = diff<0>(mat);
    // 输出：
    // 3 3 3
    // 3 3 3
    
    // 列方向差分（默认axis=1）
    auto diff_cols = diff(mat);
    // 输出：
    // 1 1
    // 1 1
    // 1 1
    return 0;
}

这个实现的优势：

• 依托Eigen表达式模板，无临时变量额外开销
• 编译期确定轴参数，性能和手写代码完全一致
• 自动适配所有Eigen矩阵类型（MatrixXd、Matrix<float, 3,3>等）

三、自定义加权平均类函数的实现思路 #

以行方向加权平均为例（每个元素与下一个元素按权重w和1-w计算），同样用模板实现：

template <typename Derived, typename Scalar = typename Derived::Scalar>
Eigen::Matrix<Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime - 1>
weighted_avg_row(const Eigen::MatrixBase<Derived>& mat, Scalar w = 0.5) {
    // 左列*w + 右列*(1-w)，取到倒数第二列
    return w * mat.leftCols(mat.cols() - 1) + (1 - w) * mat.rightCols(mat.cols() - 1);
}

如果要支持通用轴，同样用constexpr分支处理行/列方向，核心就是利用Eigen的块操作（leftCols、rightCols、topRows、bottomRows）获取相邻元素块，再做线性组合。

四、额外性能优化技巧 #

• 开启Eigen优化宏：编译时加上-DEIGEN_NO_DEBUG和-O3，关闭调试检查并启用最大优化
• 优先用固定大小矩阵：如果矩阵大小编译期已知（比如Matrix<float, 5,5>），Eigen会生成更高效的代码，避免动态内存分配
• 避免不必要拷贝：尽量返回Eigen表达式类型（或用auto推导），延迟计算减少中间变量

注意：如果需要处理3D及以上张量，Eigen原生支持有限，可考虑使用Eigen的Tensor模块，或者用模板递归处理维度，但多数有限元/差分场景下2D矩阵足够覆盖需求。

内容的提问来源于stack exchange，提问作者Ash Pera