Eigen行列をループする最も効率的な方法

Question

私は、負の数と正の数の「分けられた和」、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に存在していますか？）を使用するような他の方法がありますか？

Answer 1

Eigenは、デフォルトでカラムメジャー（Fortran）の順番で行列を割り当てます（documentation）。

行列を反復する最速の方法は、順番が間違っていると、キャッシュミスの数が増えます（行列がL1に収まらない場合、計算時間が支配的になるため、あなたの計算時間）をキャッシュライン/エレメンツ（おそらく64/8 = 8）の係数で除算します。

あなたのマトリックスがL1キャッシュに収まるなら、これは違いはありませんが、いいコンパイラがループをベクトル化できるはずです。AVXを有効にすると（光沢のある新しいコアi7で） 4回ほど。 （256ビット/ 64ビット）。

最後に、Eigenのビルトイン関数でスピードアップを期待することはありません（私は反復子があるとは思いませんが、間違いかもしれません）。同じ（非常に単純な）コード。

TLDR：繰り返しの順序を入れ替えるには、行インデックスを最も速く変更する必要があります。

Answer 2

xs（i、j）をループ内の一時変数に格納して、関数を1回だけ呼び出すようにしてください。

Answer 3

私が速くなる方法チェックアウトするためにいくつかのベンチマークをした、私は（秒）以下の結果を得た：

12 
30 
3 
6 
23 
3 

@jleahyにより示唆されるように、最初の行は、繰り返しをやっています。 2番目の行は、私のコード（@jleahyの逆順）のコードで行ったように、繰り返しを行っています。 3行目はPlainObjectBase::data()を使用してこれを反復しています。これはfor (int i = 0; i < matrixObject.size(); i++)です。 他の3行は、上記と同じですが、@ lucas92のように一時的に繰り返されます。

私は同じテストを行いましたが、/ 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; 
} 

Answer 4

このコードは、マトリックスのすべてのエントリの合計、つまり

return xs.sum(); 

私はそれが単一パスのみですので、それがより良い実行することになり引き受ける、さらに固有値は、最適なパフォーマンスのためにパスを配置する方法を「知っている」べき：、あなたはこれを行うことができます。

、しかし、あなたは二つのパスを保持したい場合、あなたはこのように、係数ごとの減速機構を使用して、これを表現することができます：正を選択するブール係数単位の選択を使用しています

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

否定的なエントリ。私はそれがハンド・ローリング・ループよりも優れているかどうかは分かりませんが、理論的にはこの方法でコード化することで、Eigen（そしてコンパイラ）はあなたの意図についてもっと知り、結果を改善する可能性があります。