矩阵乘法的 Eigen 中意外的大运行时变化

Question

我正在比较在 Eigen 中执行等效矩阵运算的方法，并且得到了非常不同的运行时间，包括一些非直观的结果。

我正在比较矩阵乘法的三种数学等效形式：

wx * 转置（数据）

我比较的三种形式是：

result = wx * data.transpose()（直接乘法版本）
result.noalias() = wx * data.transpose()（无别名版本）
result = (data * wx.transpose()).transpose()（转置版本）

我也在使用列主存储和行主存储进行测试。

使用列主存储，转置版本比直接乘法和无别名版本快得多（一个数量级），两者在运行时大致相等。

使用行主存储，noalias 和转置版本在运行时都比直接乘法快得多。

我了解到 Eigen 使用惰性求值，运算返回的直接结果通常是表达式模板，而不是中间值。我也明白 matrix * matrix 操作在它们是右侧的最后一个操作时总是会产生一个临时的，以避免别名问题，因此我试图通过 noalias() 来加快速度。

我的主要问题：

为什么转置版本总是明显更快，即使（在列主存储的情况下）当我明确声明 noalias 所以没有创建临时文件时也是如此？
为什么在使用列主存储时，运行时的（显着）差异只出现在直接乘法和 noalias 版本之间？

我为此使用的代码如下。它正在使用 gcc 4.9.2 编译，在 Centos 6 安装上，使用以下命令行。

g++ eigen_test.cpp -O3 -std=c++11 -o eigen_test -pthread -fopenmp -finline-functions

using Matrix = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor>;
  // using Matrix = Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
  int wx_rows = 8000;
  int wx_cols = 1000;
  int samples = 1;
  // Eigen::MatrixXf matrix = Eigen::MatrixXf::Random(matrix_rows, matrix_cols);
  Matrix wx = Eigen::MatrixXf::Random(wx_rows, wx_cols);
  Matrix data = Eigen::MatrixXf::Random(samples, wx_cols);

  Matrix result;

  unsigned int iterations = 10000;
  float sum = 0;
  auto before = std::chrono::high_resolution_clock::now();
  for (unsigned int ii = 0; ii < iterations; ++ii)
  {
    result = wx * data.transpose();
    sum += result(result.rows() - 1, result.cols() - 1);
  }
  auto after = std::chrono::high_resolution_clock::now();
  auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(after - before).count();
  std::cout << "original sum: " << sum << std::endl;
  std::cout << "original time (ms): " << duration << std::endl;
  std::cout << std::endl;

  sum = 0;
  before = std::chrono::high_resolution_clock::now();
  for (unsigned int ii = 0; ii < iterations; ++ii)
  {
    result.noalias() = wx * data.transpose();
    sum += result(wx_rows - 1, samples - 1);
  }
  after = std::chrono::high_resolution_clock::now();
  duration = std::chrono::duration_cast<std::chrono::milliseconds>(after - before).count();
  std::cout << "alias      sum: " << sum << std::endl;
  std::cout << "alias time (ms)     : " << duration << std::endl;
  std::cout << std::endl;

  sum = 0;
  before = std::chrono::high_resolution_clock::now();
  for (unsigned int ii = 0; ii < iterations; ++ii)
  {
    result = (data * wx.transpose()).transpose();
    sum += result(wx_rows - 1, samples - 1);
  }
  after = std::chrono::high_resolution_clock::now();
  duration = std::chrono::duration_cast<std::chrono::milliseconds>(after - before).count();
  std::cout << "new      sum: " << sum << std::endl;
  std::cout << "new time (ms)     : " << duration << std::endl;

Answer 1

一半的解释是因为，在 Eigen 的当前版本中，多线程是通过将工作拆分到结果的列块上来实现的（以及正确的-手边）。只有 1 列，不会发生多线程。在列主要案例中，这解释了案例 1 和案例 2 表现不佳的原因。另一方面，情况 3 被评估为：

column_major_tmp.noalias() = data * wx.transpose();
result = column_major_tmp.transpose();

并且由于 wx.transpose().cols() 很大，多线程是有效的。

要理解行优先的情况，您还需要知道内部矩阵乘积是针对列优先目标实现的。如果目标是行优先的，如情况 2，则乘积被转置，所以真正发生的是：

row_major_result.transpose().noalias() = data * wx.transpose();

所以我们回到案例 3，但没有临时的。

这显然是当前 Eigen 对高度不平衡矩阵大小的多线程实现的限制。理想情况下，线程应该根据手头矩阵的大小分布在行块 and/or 列块上。

顺便说一句，您还应该使用 -march=native 进行编译，让 Eigen 充分利用您的 CPU（AVX、FMA、AVX512...）。

矩阵乘法的 Eigen 中意外的大运行时变化

Unexpected and large runtime variations in Eigen for matrix multiplies

eigen