Najbardziej efektywny sposób na przechodzenie przez macierz Eigen

Question

Tworzę kilka funkcji, które wykonują takie rzeczy, jak "oddzielona suma" liczby ujemnej i dodatniej, kahan, pairwise i inne, gdzie nie ma znaczenia kolejność, w jakiej biorę elementy z matrycy, na przykład:

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

Teraz, ja, aby ten był jak najbardziej skuteczny, więc moje pytanie, byłoby lepiej pętlę przez każdą z kolumn z każdego rzędu, jak wyżej, lub wykonaj czynności odwrotne niż następujące:

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

(Przypuszczam, że to zależy od kolejność przydzielania elementów macierzy w pamięci, ale nie znalazłem tego w instrukcji Eigen).

Są również inne alternatywne sposoby, takie jak używanie iteratorów (czy istnieją w Eigen?), Które mogą być nieco szybsze?

Answer 1

Eigen przydziela domyślnie macierze w porządku kolumna-major (Fortran) (documentation).

Najszybszy sposób na iterację macierzy jest w porządku przechowywania, robienie jej w niewłaściwy sposób zwiększy liczbę braków w pamięci podręcznej (które, jeśli macierz nie mieści się w L1, zdominują twój czas obliczeniowy, więc czytaj wzrost czas obliczeń) przez współczynnik pamięci podręcznej/pamięci podręcznej (prawdopodobnie 64/8 = 8).

Jeśli macierz mieści się w pamięci podręcznej L1, nie będzie to miało znaczenia, ale dobry kompilator powinien być w stanie wektoryzować pętlę, co przy włączonym AVX (na błyszczącym nowym rdzeniu i7) może dać przyspieszenie nawet 4 razy. (256 bitów/64 bity).

Nareszcie nie oczekuj, że jakiekolwiek wbudowane funkcje Eigena spowodują przyspieszenie (nie sądzę, że w każdym razie są iteratory, ale mogę się mylić), po prostu podadzą ci ten sam (bardzo prosty) kod.

TLDR: Zamień kolejność iteracji, musisz szybko zmienić indeks wiersza.

Answer 2

Spróbuj zapisać xs (i, j) wewnątrz zmiennej tymczasowej wewnątrz pętli, aby wywołać funkcję tylko raz.

Answer 3

zrobiłem jakieś odniesienia do kasy, w jaki sposób jest szybszy, mam następujące wyniki (w sekundach):

12 
30 
3 
6 
23 
3 

Pierwsza linia robi iteracji jak sugeruje @jleahy. Druga linia wykonuje iterację, tak jak to zrobiłem w moim kodzie w pytaniu (odwrotna kolejność @jleahy). Trzecia linia wykonuje iterację za pomocą PlainObjectBase::data(), tak jak ta for (int i = 0; i < matrixObject.size(); i++). Pozostałe 3 wiersze powtarzają to samo, co powyżej, ale z tymczasowym, zgodnie z sugestią @ lucas92

Wykonałem te same testy, ale używając zamiennika/jeśli jeszcze inaczej. */Dla/else/(bez specjalnego traktowania dla rzadkich matrix) i otrzymałem następujące informacje (w sekundach):

10 
27 
3 
6 
24 
2 

Ponowne wykonanie testów dało mi podobne wyniki. Użyłem g++ 4.7.3 z -O3. Kod:

#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

Zauważam, że kod jest równoważny sumie wszystkich wpisów w macierzy, tj., Można po prostu to zrobić:

return xs.sum(); 

Przypuszczam, że byłoby lepiej, ponieważ jest to tylko jeden karnet, a ponadto Eigen powinien „wiedzieć” jak zorganizować karnety dla optymalnej wydajności.

Jeśli jednak chcesz zachować dwa karnety, można wyrazić za pomocą mechanizmów redukcji współczynnika mądry, tak:

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

który wykorzystuje logiczną współczynnik mądry wybór, aby wybrać pozytywna i negatywne wpisy. Nie wiem, czy byłaby ona lepsza od pętli ręcznych, ale teoretycznie kodowanie w ten sposób pozwala Eigenowi (i kompilatorowi) dowiedzieć się więcej o twojej intencji i ewentualnie poprawić wyniki.