11

我正在创建一些函数来执行诸如负数和正数的“分离总和”、kahan、pairwise 和其他与我从矩阵中获取元素的顺序无关的事情，例如：

template <typename T, int R, int C>
inline T sum(const Eigen::Matrix<T,R,C>& xs)
{
  T sumP(0);
  T sumN(0);
  for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
   for (size_t j = 0; j < nCols; ++j)
   {
        if (xs(i,j)>0)
          sumP += xs(i,j);
        else if (xs(i,j)<0) //ignore 0 elements: improvement for sparse matrices I think
          sumN += xs(i,j);
   }
 return sumP+sumN;
}

现在，我想让它尽可能高效，所以我的问题是，最好像上面那样遍历每一行的每一列，还是像下面这样相反：

for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
  for (size_t j = 0; j < nRows; ++j)

（我想这取决于矩阵元素在内存中分配的顺序，但我在 Eigen 的手册中找不到它）。

此外，还有其他替代方法，例如使用迭代器（它们是否存在于 Eigen 中？）可能会稍微快一些？

4

4 回答 4

18

我做了一些基准测试来检查哪种方式更快，我得到了以下结果（以秒为单位）：

12
30
3
6
23
3

第一行是按照@jleahy 的建议进行迭代。第二行是像我在问题中的代码中所做的那样进行迭代（@jleahy 的逆序）。第三行是使用PlainObjectBase::data()like this进行迭代for (int i = 0; i < matrixObject.size(); i++)。其他 3 行重复与上述相同，但由@lucas92 建议临时

我也做了同样的测试，但是用 /if else.*/ 代替 /else/ （对稀疏矩阵没有特殊处理），我得到了以下结果（以秒为单位）：

10
27
3
6
24
2

再次进行测试给了我非常相似的结果。我g++ 4.7.3与-O3. 编码：

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

using namespace std;

 template <typename T, int R, int C>
    inline T sum_kahan1(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
      for (size_t j = 0; j < nRows; ++j)
      {
          if (xs(j,i)>0)
          {
            yP = xs(j,i) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if (xs(j,i)<0)
          {
            yN = xs(j,i) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan2(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
      for (size_t j = 0; j < nCols; ++j)
      {
          if (xs(i,j)>0)
          {
            yP = xs(i,j) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if (xs(i,j)<0)
          {
            yN = xs(i,j) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }


 template <typename T, int R, int C>
    inline T sum_kahan3(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
    for (size_t i = 0, size = xs.size(); i < size; i++)
      {
          if ((*(xs.data() + i))>0)
          {
            yP = (*(xs.data() + i)) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if ((*(xs.data() + i))<0)
          {
            yN = (*(xs.data() + i)) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan1t(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
      for (size_t j = 0; j < nRows; ++j)
      {
      T temporary = xs(j,i);
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if (temporary<0)
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan2t(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
      for (size_t j = 0; j < nCols; ++j)
      {
      T temporary = xs(i,j);
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if (temporary<0)
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }


 template <typename T, int R, int C>
    inline T sum_kahan3t(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
    for (size_t i = 0, size = xs.size(); i < size; i++)
      {
      T temporary = (*(xs.data() + i));
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else if (temporary<0)
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan1e(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
      for (size_t j = 0; j < nRows; ++j)
      {
          if (xs(j,i)>0)
          {
            yP = xs(j,i) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = xs(j,i) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan2e(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
      for (size_t j = 0; j < nCols; ++j)
      {
          if (xs(i,j)>0)
          {
            yP = xs(i,j) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = xs(i,j) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }


 template <typename T, int R, int C>
    inline T sum_kahan3e(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
    for (size_t i = 0, size = xs.size(); i < size; i++)
      {
          if ((*(xs.data() + i))>0)
          {
            yP = (*(xs.data() + i)) - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = (*(xs.data() + i)) - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan1te(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nCols; ++i)
      for (size_t j = 0; j < nRows; ++j)
      {
      T temporary = xs(j,i);
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }

 template <typename T, int R, int C>
    inline T sum_kahan2te(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
      for (size_t i = 0, nRows = xs.rows(), nCols = xs.cols(); i < nRows; ++i)
      for (size_t j = 0; j < nCols; ++j)
      {
      T temporary = xs(i,j);
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }


 template <typename T, int R, int C>
    inline T sum_kahan3te(const Eigen::Matrix<T,R,C>& xs) {
      if (xs.size() == 0) return 0;
      T sumP(0);
      T sumN(0);
      T tP(0);
      T tN(0);
      T cP(0);
      T cN(0);
      T yP(0);
      T yN(0);
    for (size_t i = 0, size = xs.size(); i < size; i++)
      {
      T temporary = (*(xs.data() + i));
          if (temporary>0)
          {
            yP = temporary - cP;
          tP = sumP + yP;
          cP = (tP - sumP) - yP;
          sumP = tP;
          }
        else
          {
            yN = temporary - cN;
          tN = sumN + yN;
          cN = (tN - sumN) - yN;
          sumN = tN;
          }
      }
      return sumP+sumN;
    }


int main() {

    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> test = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>::Random(10000,10000);

    cout << "start" << endl;   
    int now;

    now = time(0);
    sum_kahan1(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan2(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan3(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan1t(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan2t(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan3t(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan1e(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan2e(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan3e(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan1te(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan2te(test);
    cout << time(0) - now << endl;

    now = time(0);
    sum_kahan3te(test);
    cout << time(0) - now << endl;

    return 0;
}
于 2013-04-29T20:50:23.923 回答
11

Eigen 默认以列优先 (Fortran) 顺序分配矩阵（文档）。

遍历矩阵的最快方法是按存储顺序，以错误的方式执行会增加缓存未命中的数量（如果您的矩阵不适合 L1 将主导您的计算时间，因此读取会增加您的计算时间）缓存线/elemsize 的一个因子（可能是 64/8=8）。

如果您的矩阵适合 L1 缓存，这不会有什么不同，但一个好的编译器应该能够对循环进行矢量化，启用 AVX（在闪亮的新核心 i7 上）可以让您的加速高达 4 倍. （256 位/64 位）。

最后，不要指望 Eigen 的任何内置函数都能加快速度（无论如何我认为没有迭代器，但我可能弄错了），它们只会给你同样的效果（非常简单） 代码。

TLDR：交换你的迭代顺序，你需要最快地改变行索引。

于 2013-04-29T19:22:24.717 回答
4

我注意到代码相当于矩阵中所有条目的总和，即，您可以这样做：

return xs.sum();

我认为这会表现得更好，因为它只是一次通过，此外，Eigen 应该“知道”如何安排通过以获得最佳性能。

但是，如果您想保留这两个通道，则可以使用系数缩减机制来表达这一点，如下所示：

return (xs.array() > 0).select(xs, 0).sum() +
       (xs.array() < 0).select(xs, 0).sum();

它使用布尔系数选择来选择正负条目。我不知道它是否会胜过手动循环，但理论上以这种方式编码可以让 Eigen（和编译器）更多地了解您的意图，并可能改善结果。

于 2015-08-20T01:10:48.283 回答
1

尝试将 xs(i,j) 存储在循环内的临时变量中，以便您只调用该函数一次。

于 2013-04-29T17:15:46.860 回答