Eigen：编码风格对性能的影响

Question

据我所了解的Eigen（这里）,似乎operator=()作为一种“障碍”，可以防止惰性评估——例如，它会导致Eigen停止返回表达式模板，并实际执行（优化）计算，将结果存储到=的左侧。

这似乎意味着一个人的“编码风格”对性能有影响——即使用命名变量来存储中间计算结果可能会对性能产生负面影响，因为这会导致某些部分的计算“过早地”被评估。

为了验证我的直觉，我写了一个示例，并对结果感到惊讶（完整代码请看这里）：

using ArrayXf  = Eigen::Array <float, Eigen::Dynamic, Eigen::Dynamic>;
using ArrayXcf = Eigen::Array <std::complex<float>, Eigen::Dynamic, Eigen::Dynamic>;

float test1( const MatrixXcf & mat )
{
    ArrayXcf arr  = mat.array();
    ArrayXcf conj = arr.conjugate();
    ArrayXcf magc = arr * conj;
    ArrayXf  mag  = magc.real();
    return mag.sum();
}

float test2( const MatrixXcf & mat )
{
    return ( mat.array() * mat.array().conjugate() ).real().sum();
}

float test3( const MatrixXcf & mat )
{
    ArrayXcf magc   = ( mat.array() * mat.array().conjugate() );

    ArrayXf mag     = magc.real();
    return mag.sum();
}

上述内容提供了三种计算复数矩阵系数幅值的不同方法。

1. `test1` 可以说是一步步地进行每个部分的计算。 2. `test2` 则通过一个表达式完成整个计算。 3. `test3` 采用了“混合”方法——使用了一定量的中间变量。

我本来期望由于 `test2` 将整个计算打包到一个表达式中，Eigen 应该能够利用它并对整个计算进行全局优化，从而提供最佳性能。

然而，结果却出人意料（显示的数字是每个测试执行1000次的总微秒数）。

test1_us: 154994
test2_us: 365231
test3_us: 36613

（这是使用g++ -O3编译的——详见gist。）

我原以为最快的版本（test2）会是最慢的。 而我预计最慢的版本（test1）居然排在中间。

因此，我的问题是：

1. 为什么test3比其他方案表现得更好？
2. 除了深入研究汇编代码之外，有没有一种技术可以让人们了解Eigen实际上如何实现您的计算？
3. 有没有一组准则可以遵循，以在您的Eigen代码中平衡性能和可读性（使用中间变量）？

在更复杂的计算中，将所有内容都写在一个表达式中可能会影响可读性，因此我希望找到正确的编写既易读又高效的代码的方法。

- jeremytrimble

---

我不是一个优化专家，但我会怀疑你的结果，因为你使用了-O3编译选项，并且没有捕获计算结果。完全有可能优化器会认识到funcN()没有副作用并且优化掉整个计算过程。我相信您可以使用volatile来辅助微基准测试。相关的Stack Overflow问题链接 - Steve Lorimer

请注意，使用最新的编译器时，程序总是会中止。它只能通过旧版本的编译器，因为调用的 abs 版本是整数版本... - Marc Glisse

Answer 1

看起来是GCC的一个问题。Intel编译器会给出期望的结果。

$ g++ -I ~/program/include/eigen3 -std=c++11 -O3 a.cpp -o a && ./a
test1_us: 200087
test2_us: 320033
test3_us: 44539

$ icpc -I ~/program/include/eigen3 -std=c++11 -O3 a.cpp -o a && ./a
test1_us: 214537
test2_us: 23022
test3_us: 42099

与 icpc 版本相比，gcc 似乎存在优化 test2 的问题。

为了得到更精确的结果，您可以通过执行 -DNDEBUG 关闭调试断言，如此处所示：here。

编辑：

问题1：

@ggael 给出了一个很好的答案，即 gcc 无法将求和循环向量化。我的实验也发现，对于 gcc 和 icc，test2 的速度与手写的 naive for-loop 相同，这表明向量化是原因，并且使用下面提到的方法未检测到 test2 中的临时内存分配，说明 Eigen 正确地评估了表达式。

问题2：

避免中间内存是 Eigen 使用表达式模板的主要目的。因此，Eigen 提供了一个宏 EIGEN_RUNTIME_NO_MALLOC 和一个简单的函数来使您能够在计算表达式时检查是否分配了中间内存。您可以在这里找到示例代码。请注意，这可能仅在调试模式下起作用。

EIGEN_RUNTIME_NO_MALLOC - 如果定义了，将引入一个新的开关，可以通过调用 set_is_malloc_allowed(bool) 打开和关闭。如果不允许分配内存，并且 Eigen 试图动态分配内存，会导致断言失败。默认情况下未定义。

问题3：

有一种方法可以使用中间变量并同时获得惰性求值/表达式模板引入的性能改进。

该方法是使用正确数据类型的中间变量。您应该使用表达式类型 Eigen::MatrixBase/ArrayBase/DenseBase，而不是 Eigen::Matrix/Array，以指示表达式仅被缓冲而不被评估。这意味着您应该将表达式存储为中间结果，而不是表达式的结果，前提是该中间结果在以下代码中只使用一次。

由于确定表达式类型 Eigen::MatrixBase/... 中的模板参数可能很繁琐，因此您可以使用 auto 代替。您可以在此页面上找到何时应该/不应该使用auto/表达式类型的一些提示。另一页也告诉您如何将表达式作为函数参数传递而不进行评估。

根据 @ggael 在答案中关于 .abs2() 的示例实验，我认为另一个指导原则是避免重复造轮子。

- kangshiyin

---

auto 应该也缓存中间结果。 - Yakk - Adam Nevraumont

它们应该仅用于缓冲中间值，这些中间值您不希望被计算。 - kangshiyin

Answer 2

由于进行了 .real() 步骤，Eigen 不会显式地对 test2 进行向量化。因此，它将调用标准的 complex::operator* 操作符，但不幸的是，这个操作符在 gcc 中永远不会被内联。另一方面，其他版本则使用了 Eigen 自己向量化实现的复数乘积。

相比之下，ICC 在内联复数运算符 complex::operator* 方面做得更好，因此使得 test2 在 ICC 上运行最快。您也可以将 test2 重写为：

return mat.array().abs2().sum();

为了在所有编译器上获得更好的性能：

gcc:
test1_us: 66016
test2_us: 26654
test3_us: 34814

icpc:
test1_us: 87225
test2_us: 8274
test3_us: 44598

clang:
test1_us: 87543
test2_us: 26891
test3_us: 44617

ICC在这种情况下得分非常高，这归功于其聪明的自动矢量化引擎。

解决gcc内联失败的另一种方法是，在不修改test2的情况下为complex<float>定义自己的operator*。例如，在文件顶部添加以下内容：

namespace std {
  complex<float> operator*(const complex<float> &a, const complex<float> &b) {
    return complex<float>(real(a)*real(b) - imag(a)*imag(b), imag(a)*real(b) + real(a)*imag(b));
  }
}

然后我得到：

gcc:
test1_us: 69352
test2_us: 28171
test3_us: 36501

icpc:
test1_us: 93810
test2_us: 11350
test3_us: 51007

clang:
test1_us: 83138
test2_us: 26206
test3_us: 45224

当然，这种技巧并不总是建议使用，因为与流畅版本相比，它可能会导致溢出或数值取消问题，但这正是icpc和其他向量化版本所计算的。

- ggael

---

好的，看起来我的直觉是正确的，除了gcc的一个怪癖。谢谢你的回复。至于我另外两个问题，你有什么技巧建议可以深入了解Eigen选择优化表达式的方式吗？此外，有没有办法将计算分成多个子表达式而不会妨碍惰性求值（例如，我是否可以将arr、conj、magc等声明为某种不同类型，以便稍后进行评估）？ - jeremytrimble

为什么.real()不能向量化，这是库中的缺陷吗？ - Yakk - Adam Nevraumont

使用 -fcx-limited-range（包括在 -ffast-math 中）或 -fcx-fortran-rules，gcc 将内联复数乘法。Icc 默认为不安全模式，这是一个值得质疑的选择... - Marc Glisse

实际上，使用 -fcx-... 选项中的任何一个，gcc 给出的结果与 icc 相同。 - Marc Glisse

@Yakk: 对于MatrixXd::real()的向量化很容易，但是当它在一个复杂表达式中被调用时，这就复杂多了，因为您必须请求两个数据包才能输出一个... 您因此增加了寄存器压力，并且还防止编译器删除无用代码（导致虚部）。更好的策略是将此信息传播到子表达式，但以合理的编译时间执行这个操作并不容易。 - ggael

@MarcGlisse，说得好。-fcx-*选项确实比我的重载技巧更干净。不幸的是，它们不受clang支持，而-ffast-math无法帮助clang内联复杂乘法。 - ggael

Answer 3

我曾经做过的一件事情就是经常使用auto关键字。需要记住的是，大多数Eigen表达式返回特殊的表达式数据类型（例如CwiseBinaryOp），将其赋值回一个Matrix可能会强制求值表达式（这就是你看到的）。使用auto可以让编译器推断返回类型为任何表达式类型，尽可能避免求值：

float test1( const MatrixXcf & mat )
{
    auto arr  = mat.array();
    auto conj = arr.conjugate();
    auto magc = arr * conj;
    auto mag  = magc.real();
    return mag.sum();
}

这应该更接近于你的第二个测试案例。在某些情况下，我在保持可读性的同时获得了良好的性能提升（您不需要拼写表达式模板类型）。当然，具体情况因人而异，因此请仔细进行基准测试 :)

- mindriot

Answer 4

我只是想提醒你，你的性能分析方法并不是最优的，所以实际上问题可能只是你的性能分析方法。

由于有许多像缓存局部性这样需要考虑的因素，你应该按照以下方式进行性能分析：

int warmUpCycles = 100;
int profileCycles = 1000;

// TEST 1
for(int i=0; i<warmUpCycles ; i++)
      doTest1();

auto tick = std::chrono::steady_clock::now();
for(int i=0; i<profileCycles ; i++)
      doTest1();  
auto tock = std::chrono::steady_clock::now();
test1_us = (std::chrono::duration_cast<std::chrono::microseconds>(tock-tick)).count(); 

// TEST 2


// TEST 3

一旦你以正确的方式进行了测试，那么你就可以得出结论。

我非常怀疑，由于您一次只对一个操作进行分析，因此在第三次测试时会使用缓存版本，因为编译器很可能重新排序操作。

此外，您应该尝试不同的编译器，以查看问题是否是模板展开（模板优化有深度限制：很可能您可以用单个大表达式击中它）。

此外，如果Eigen支持移动语义，那么没有理由认为一个版本应该更快，因为并不总是保证可以优化表达式。

请尝试并告诉我，这很有趣。还要确保使用像-O3这样的标志启用了优化，没有优化的分析是没有意义的。

为了防止编译器将所有内容都优化掉，请使用来自文件或cin的初始输入，然后重新将输入馈送到函数中。

- CoffeDeveloper

---

移动语义在这里没有帮助，因为临时对象会引入缓存未命中，并阻止通过融合操作实现的优化机会。 - ggael