Eigen C++：稀疏矩阵操作的性能答案

【问题标题】：Eigen C++: Performance of sparse-matrix manipulationsEigen C++：稀疏矩阵操作的性能
【发布时间】：2016-08-03 21:35:01
【问题描述】：

谁能解释特征稀疏矩阵的以下行为？我一直在研究aliasing 和lazy evaluation，但我似乎无法改进这个问题。技术规格：我在 Ubuntu 16.10 上使用最新的 Eigen 稳定版本，带有 g++ 编译器，没有优化标志。

假设我用以下方式定义一个简单的身份：

SparseMatrix<double> spIdent(N,N);
spIdent.reserve(N);
spIdent.setIdentity();

然后用它执行这些操作

spIdent-spIdent;
spIdent*spIdent;
spIdent - spIdent*spIdent;

并测量所有三个的计算时间。我得到的是这个

0 Computation time: 2.6e-05
1 Computation time: 2e-06 
2 Computation time: 1.10706

这意味着任何一个操作都很快，但组合起来非常慢。 noalias() 方法只为 dense 矩阵定义，而且在我的密集示例中，它并没有太大的区别。有什么启示吗？

MCVE：

#include <iostream>
#include <ctime>
#include "../Eigen/Sparse"

using namespace std;
using namespace Eigen;

int main() {

unsigned int N=2000000;

SparseMatrix<double> spIdent(N,N);
spIdent.reserve(N);
spIdent.setIdentity();

clock_t start=clock();
spIdent*spIdent;
cout << "0 Computation time: " << float(clock() - start)/1e6 << '\n';

start=clock();
spIdent-spIdent;
cout << "1 Computation time: " << float(clock() - start)/1e6 << '\n';

start=clock();
spIdent - (spIdent*spIdent);
cout << "2 Computation time: " << float(clock() - start)/1e6 << '\n';

return 0;

}

【问题讨论】：

呃，使用优化标志？
问题不只是让它更快......而是理解为什么第三种情况如此缓慢！
前两个可能优化为根本没有代码。第三个不是。
在没有开启优化的情况下询问速度，是徒劳的。
你能创建一个 MCVE 代码来重现你的计时结果吗？ stackoverflow.com/help/mcve

标签： c++ matrix eigen

【解决方案1】：

它并没有像惰性求值那样得到优化，嗯，非常惰性。看产品。被调用的代码是（至少在这台机器上包含的任何版本的 Eigen 中）：

template<typename Derived>
template<typename OtherDerived>
inline const typename SparseSparseProductReturnType<Derived,OtherDerived>::Type
SparseMatrixBase<Derived>::operator*(const SparseMatrixBase<OtherDerived> &other) const
{
  return typename SparseSparseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
}

返回产品的表达式（即惰性）。这个表达式没有做任何事情，所以成本为零。差异也是如此。现在，在执行a-a*a 时，a*a 是一个表达式。然后它遇到operator-。这看到了右手边的表达式。然后将此表达式评估为临时（即花费时间）以便在operator- 中使用它。为什么评价为临时？阅读this 了解他们的逻辑（搜索“第二种情况”）。

operator- 是一个CwiseBinaryOp，右侧是产品表达式。 CwiseBinaryOp 做的第一件事是将右手边分配给成员：

EIGEN_STRONG_INLINE CwiseBinaryOp(const Lhs& aLhs, const Rhs& aRhs, const BinaryOp& func = BinaryOp())
      : m_lhs(aLhs), m_rhs(aRhs), m_functor(func)

(m_rhs(aRhs)) 依次调用SparseMatrix 构造函数：

/** Constructs a sparse matrix from the sparse expression \a other */
template<typename OtherDerived>
inline SparseMatrix(const SparseMatrixBase<OtherDerived>& other)
  : m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0)
{
  ...
  *this = other.derived();
}

这反过来又调用operator=，这将（如果我错了，有人纠正我）总是触发评估，在这种情况下，是临时的。

【讨论】：

那么，是否可以避免表达式a - a * a中的临时变量？
使用 Eigen 3.2 中采用的自下而上的方法，我不这么认为。 Eigen 3.3（测试版）可能采用自上而下的方法（我想我记得在某处看到过，但我可能是错的）。如果你问理论上是否可行，那么当然可以。
更准确地说，在 Eigen 3.2 中，第三条语句中的产品立即被评估为临时产品。 Eigen 3.3 不再是这种情况。因此，对于 Eigen 3.3，所有三个语句都归结为无操作。但是，即使使用 3.3，一旦将第三条语句分配给给定的矩阵，那么产品将在临时范围内进行评估。通过点积懒惰地评估它会非常慢。在这种情况下，临时的是值得的！
感谢您的澄清。

【解决方案2】：

好吧，正如人们所提到的，在前两个语句中，代码被完全优化掉了（我已经用当前版本的 g++ 和-O3 set 进行了测试）。反汇编显示了第二个语句：

  400e78:   e8 03 fe ff ff          callq  400c80 <clock@plt>   # timing begins
  400e7d:   48 89 c5                mov    %rax,%rbp
  400e80:   e8 fb fd ff ff          callq  400c80 <clock@plt>   # timing ends

Fop 第三部分实际上发生了一些事情，调用了 Eigen 库代码：

  400ede:   e8 9d fd ff ff          callq  400c80 <clock@plt>   # timing begins
  400ee3:   48 89 c5                mov    %rax,%rbp
  400ee6:   8b 44 24 58             mov    0x58(%rsp),%eax
  400eea:   39 44 24 54             cmp    %eax,0x54(%rsp)
  400eee:   c6 44 24 20 00          movb   $0x0,0x20(%rsp)
  400ef3:   48 89 5c 24 28          mov    %rbx,0x28(%rsp)
  400ef8:   48 89 5c 24 30          mov    %rbx,0x30(%rsp)
  400efd:   48 c7 44 24 38 00 00    movq   $0x0,0x38(%rsp)
  400f04:   00 00 
  400f06:   c6 44 24 40 01          movb   $0x1,0x40(%rsp)
  400f0b:   0f 85 99 00 00 00       jne    400faa <main+0x22a>
  400f11:   48 8d 4c 24 1f          lea    0x1f(%rsp),%rcx
  400f16:   48 8d 54 24 20          lea    0x20(%rsp),%rdx
  400f1b:   48 8d bc 24 90 00 00    lea    0x90(%rsp),%rdi
  400f22:   00 
  400f23:   48 89 de                mov    %rbx,%rsi
  400f26:   e8 25 1a 00 00          callq  402950 <_ZN5Eigen13CwiseBinaryOpINS_8internal20scalar_difference_opIdEEKNS_12SparseMatrixIdLi0EiEEKNS_19SparseSparseProductIRS6_S8_EEEC1ES8_RSA_RKS3_>
  400f2b:   48 8d bc 24 a0 00 00    lea    0xa0(%rsp),%rdi
  400f32:   00 
  400f33:   e8 18 02 00 00          callq  401150 <_ZN5Eigen12SparseMatrixIdLi0EiED1Ev>
  400f38:   e8 43 fd ff ff          callq  400c80 <clock@plt>   # timing ends

我猜在这种情况下编译器无法确定计算结果没有被使用，与前两种情况相反。

如果您查看documentation，您会发现像+ 这样对稀疏矩阵的简单操作返回的不是矩阵，而是代表结果的CwiseUnaryOp。我想如果你不在某个地方使用这个类，那么生成的矩阵永远不会被构造出来。

【讨论】：

【解决方案3】：

我认为正如@hfhc2 所提到的，代码中的前两个语句完全由编译器优化（因为其余的不需要结果）。在第三个语句中，最有可能产生一个辅助中间变量来存储spIdent*spIdent 的临时结果。为了清楚地看到这一点，请考虑以下包含显式复制分配的示例：

#include <iostream>
#include <ctime>
#include <Eigen/Sparse>

using namespace std;
using namespace Eigen;

int main () {

   const unsigned int N = 2000000;

   SparseMatrix<double> spIdent(N,N);
   SparseMatrix<double> a(N,N), b(N,N), c(N,N);

   spIdent.reserve(N);
   spIdent.setIdentity();

   clock_t start = clock();
   a = spIdent*spIdent;
   cout << "0 Computation time: " << float(clock() - start)/1e6 << endl;

   start = clock();
   b = spIdent-spIdent;
   cout << "1 Computation time: " << float(clock() - start)/1e6 << endl;

   start = clock();
   c = a - b;
   cout << "2 Computation time: " << float(clock() - start)/1e6 << endl;

   return 0;

}

测量的时间（没有编译器优化）是 [for openSUSE 12.2 (x86_64), g++ 4.7.1, Intel 2core 2GHz CPU]：

0 Computation time: 1.58737
1 Computation time: 0.417798
2 Computation time: 0.428174

这似乎很合理。

【讨论】：