Cách hiệu quả nhất để lặp qua ma trận Eigen

Question

Tôi đang tạo một số chức năng để thực hiện những việc như "số tiền được tách biệt" của số âm và dương, kahan, pairwise và các thứ khác, không quan trọng thứ tự tôi lấy các yếu tố từ các ma trận, ví dụ:

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

Bây giờ, tôi muốn thực hiện điều này như hiệu quả càng tốt, vì vậy câu hỏi của tôi là, nó sẽ là tốt hơn để lặp qua mỗi cột của mỗi dòng như trên, hoặc làm ngược lại như sau:

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

(Tôi cho rằng điều này phụ thuộc vào thứ tự các phần tử ma trận được cấp phát trong bộ nhớ, nhưng tôi không thể tìm thấy điều này trong sách hướng dẫn của Eigen).

Ngoài ra, có các cách thay thế khác như sử dụng trình vòng lặp (chúng có tồn tại trong Eigen không?) Có thể nhanh hơn một chút?

Answer 1

Eigen phân bổ ma trận theo thứ tự cột lớn (Fortran) theo mặc định (documentation). Cách nhanh nhất để lặp qua ma trận là theo thứ tự lưu trữ, làm theo cách sai vòng sẽ làm tăng số lần nhớ cache (nếu ma trận của bạn không phù hợp với L1 sẽ thống trị thời gian tính toán của bạn, vì vậy hãy đọc tăng lên) thời gian tính toán của bạn) bởi một yếu tố của cacheline/elemsize (có thể là 64/8 = 8).

Nếu ma trận của bạn phù hợp với bộ đệm L1, điều này sẽ không tạo ra sự khác biệt, nhưng trình biên dịch tốt sẽ có thể vectơ vòng lặp, với AVX được bật (trên một lõi mới sáng bóng i7) có thể giúp bạn tăng tốc nhiều gấp 4 lần. (256 bit/64 bit). Cuối cùng, tôi không mong đợi bất kỳ chức năng tích hợp nào của Eigen để giúp bạn tăng tốc (tôi không nghĩ có lặp lại dù sao đi nữa, nhưng tôi có thể nhầm lẫn), họ sẽ chỉ cung cấp cho bạn cùng một mã (rất đơn giản).

TLDR: Hoán đổi thứ tự lặp lại của bạn, bạn cần phải thay đổi chỉ mục hàng nhanh nhất.

Answer 2

Cố gắng lưu trữ xs (i, j) bên trong một biến tạm thời bên trong vòng lặp để bạn chỉ gọi hàm này một lần.

Answer 3

tôi đã làm một số tiêu chuẩn để kiểm tra đó là cách nhanh hơn, tôi nhận được kết quả như sau (tính bằng giây):

12 
30 
3 
6 
23 
3 

Dòng đầu tiên được thực hiện lặp đi lặp lại theo đề nghị của @jleahy. Dòng thứ hai đang thực hiện lặp lại như tôi đã thực hiện trong mã của tôi trong câu hỏi (thứ tự ngược của @jleahy). Dòng thứ ba đang thực hiện lặp lại sử dụng PlainObjectBase::data() như thế này for (int i = 0; i < matrixObject.size(); i++). 3 dòng còn lại lặp lại giống như trên, nhưng tạm thời được gợi ý bởi @ lucas92

Tôi cũng đã làm các thử nghiệm tương tự nhưng sử dụng thay thế/nếu khác. */For/else/(không xử lý đặc biệt cho thưa thớt ma trận) và tôi nhận được thông tin sau (tính bằng giây):

10 
27 
3 
6 
24 
2 

Thực hiện các kiểm tra lại cho kết quả tương tự. Tôi đã sử dụng g++ 4.7.3 với -O3. Mã:

#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

Tôi nhận thấy rằng mã này tương đương với tổng của tất cả các mục nhập trong ma trận, tức là, Bạn chỉ có thể làm điều này:

return xs.sum(); 

Tôi cho rằng sẽ thực hiện tốt hơn kể từ khi nó chỉ là một pass duy nhất, và hơn nữa Eigen phải "biết" làm thế nào để sắp xếp các đường chuyền cho hiệu suất tối ưu.

Tuy nhiên, nếu bạn muốn giữ lại hai đèo, bạn có thể thể hiện điều này bằng cách sử dụng cơ chế giảm hệ số khôn ngoan, như thế này:

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

trong đó sử dụng boolean hệ số khôn ngoan chọn để chọn tích cực và các mục tiêu cực. Tôi không biết liệu nó có hoạt động tốt hơn các vòng quay tay hay không, nhưng theo lý thuyết mã hóa nó theo cách này cho phép Eigen (và trình biên dịch) biết nhiều hơn về ý định của bạn, và có thể cải thiện kết quả.