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

Question

从我读到的关于 Eigen (here) 的内容来看，operator=() 似乎充当了懒惰评估的“障碍”——例如它会导致 Eigen 停止返回表达式模板并实际执行（优化的）计算，将结果存储到 = 的左侧。

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

为了验证我的直觉，我写了一个例子，结果让我大吃一惊 (full code here)：

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();
}

上面给出了 3 种不同的方法来计算复值矩阵中的系数大小之和。

1. test1 有点像“一次一步”进行计算的每个部分。
2. test2 在一个表达式中完成整个计算。
3. test3 采用“混合”方法——使用一些中间变量。

我有点期待，由于 test2 将整个计算打包到一个表达式中，Eigen 将能够利用这一点并全局优化整个计算，从而提供最佳性能。

但是，结果令人惊讶（显示的数字是每个测试执行 1000 次的总微秒数）：

test1_us: 154994
test2_us: 365231
test3_us: 36613

（这是使用 g++ -O3 编译的——有关详细信息，请参阅 the gist。）

我预计最快的版本 (test2) 实际上是最慢的。另外，我预计最慢的版本（test1）实际上在中间。

所以，我的问题是：

1. 为什么test3 的性能比替代品好很多？
2. 是否可以使用一种技术（无需深入研究汇编代码）来了解 Eigen 是如何实际实现您的计算的？
3. 是否有一套指导方针可以在您的 Eigen 代码中实现性能和可读性（使用中间变量）之间的良好权衡？

在更复杂的计算中，在一个表达式中执行所有操作可能会妨碍可读性，因此我有兴趣找到编写可读性和高性能代码的正确方法。

Answer 1

看起来像是 GCC 的问题。英特尔编译器给出了预期的结果。

$ 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 无法对 sum 循环进行矢量化。我的实验还发现test2和手写的naive for-loop一样快，gcc和icc都有，提示向量化是原因，在test2中没有检测到临时内存分配下面提到的方法，表明 Eigen 正确地评估了表达式。

关于问题 2

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

EIGEN_RUNTIME_NO_MALLOC - 如果已定义，则会引入一个新开关，该开关 可以通过调用 set_is_malloc_allowed(bool) 打开和关闭。如果 不允许 malloc 并且 Eigen 尝试动态分配内存 无论如何，断言失败的结果。默认情况下未定义。

关于问题 3

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

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

由于在表达式类型Eigen::MatrixBase/... 中确定模板参数可能会很痛苦，因此您可以改用auto。您可以在 this page 中找到关于何时应该/不应该使用 auto/表达式类型的一些提示。 Another page 还告诉您如何将表达式作为函数参数传递而不计算它们。

根据@ggael 的回答中关于.abs2() 的指导性实验，我认为另一个准则是避免重新发明轮子。

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 在这种情况下的优异成绩是由于其巧妙的自动矢量化引擎。

另一种在不修改test2 的情况下解决gcc 内联失败的方法是为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

当然，并不总是推荐使用此技巧，因为与 glib 版本相比，它可能会导致溢出或数值取消问题，但无论如何 icpc 和其他矢量化版本都会计算这一点。

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();
}

这应该更接近您的第二个测试用例。在某些情况下，我在保持可读性的同时获得了很好的性能改进（您确实不希望必须拼出表达式模板类型）。当然，您的里程可能会有所不同，因此请仔细进行基准测试:)

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 的初始输入，然后在函数内重新输入输入。