2016-03-28 317 views
3

我正在寫一個庫,我想要有一些基本的NxN矩陣功能,沒有任何依賴關係,它是一個學習項目。我將我的表現與Eigen進行比較。我已經能夠相當平等,甚至可以在SSE2的前面擊敗它的表現,並且AVX2在很多方面都打敗了它(它只使用SSE2,所以不會超級意外)。計算行列式的最快方法?

我的問題是我使用高斯消去來創建一個上對角化矩陣,然後乘以對角線得到行列式。我擊敗Eigen爲N < 300,但之後,Eigen吹走了我,並作爲矩陣變得更糟變得更大。鑑於所有的內存是按順序訪問的,編譯器的反彙編看起來並不可怕,我不認爲這是一個優化問題。還有更多優化可以完成,但時序看起來更像是算法時序複雜性問題,或者我沒有看到主要的SSE優勢。簡單地展開這些循環對於我來說沒有太大的幫助。

是否有更好的算法來計算行列式?

標量代碼

/* 
    Warning: Creates Temporaries! 
*/ 
template<typename T, int ROW, int COLUMN> MML_INLINE T matrix<T, ROW, COLUMN>::determinant(void) const 
{ 
    /* 
    This method assumes square matrix 
    */ 
    assert(row() == col()); 
    /* 
    We need to create a temporary 
    */ 
    matrix<T, ROW, COLUMN> temp(*this); 
    /*We convert the temporary to upper triangular form*/ 
    uint N = row(); 
    T det = T(1); 
    for (uint c = 0; c < N; ++c) 
    { 
     det = det*temp(c,c); 
     for (uint r = c + 1; r < N; ++r) 
     { 
      T ratio = temp(r, c)/temp(c, c); 
      for (uint k = c; k < N; k++) 
      { 
       temp(r, k) = temp(r, k) - ratio * temp(c, k); 
      } 
     } 
    } 

    return det; 
} 

AVX2

template<> float matrix<float>::determinant(void) const 
{ 
    /* 
    This method assumes square matrix 
    */ 
    assert(row() == col()); 
    /* 
    We need to create a temporary 
    */ 
    matrix<float> temp(*this); 
    /*We convert the temporary to upper triangular form*/ 
    float det = 1.0f; 

    const uint N = row(); 
    const uint Nm8 = N - 8; 
    const uint Nm4 = N - 4; 

    uint c = 0; 
    for (; c < Nm8; ++c) 
    { 
     det *= temp(c, c); 
     float8 Diagonal = _mm256_set1_ps(temp(c, c)); 

     for (uint r = c + 1; r < N;++r) 
     { 
      float8 ratio1 = _mm256_div_ps(_mm256_set1_ps(temp(r,c)), Diagonal); 

      uint k = c + 1; 
      for (; k < Nm8; k += 8) 
      { 
       float8 ref = _mm256_loadu_ps(temp._v + c*N + k); 
       float8 r0 = _mm256_loadu_ps(temp._v + r*N + k); 

       _mm256_storeu_ps(temp._v + r*N + k, _mm256_fmsub_ps(ratio1, ref, r0)); 
      } 

      /*We go Scalar for the last few elements to handle non-multiples of 8*/ 
      for (; k < N; ++k) 
      { 
       _mm_store_ss(temp._v + index(r, k), _mm_sub_ss(_mm_set_ss(temp(r, k)), _mm_mul_ss(_mm256_castps256_ps128(ratio1),_mm_set_ss(temp(c, k))))); 
      } 
     } 
    } 

    for (; c < Nm4; ++c) 
    { 
     det *= temp(c, c); 
     float4 Diagonal = _mm_set1_ps(temp(c, c)); 

     for (uint r = c + 1; r < N; ++r) 
     { 
      float4 ratio = _mm_div_ps(_mm_set1_ps(temp[r*N + c]), Diagonal); 
      uint k = c + 1; 
      for (; k < Nm4; k += 4) 
      { 
       float4 ref = _mm_loadu_ps(temp._v + c*N + k); 
       float4 r0 = _mm_loadu_ps(temp._v + r*N + k); 

       _mm_storeu_ps(temp._v + r*N + k, _mm_sub_ps(r0, _mm_mul_ps(ref, ratio))); 
      } 

      float fratio = _mm_cvtss_f32(ratio); 
      for (; k < N; ++k) 
      { 
       temp(r, k) = temp(r, k) - fratio*temp(c, k); 
      } 
     } 
    } 

    for (; c < N; ++c) 
    { 
     det *= temp(c, c); 
     float Diagonal = temp(c, c); 
     for (uint r = c + 1; r < N; ++r) 
     { 
      float ratio = temp[r*N + c]/Diagonal; 
      for (uint k = c+1; k < N;++k) 
      { 
       temp(r, k) = temp(r, k) - ratio*temp(c, k); 
      } 
     } 
    } 

    return det; 
} 
+1

鑑於您已經準備好深入細節,並且Eigen是開源的... [爲什麼不看它的功能](https://github.com/RLovelett/eigen/search?utf8 =%E2%9C%93&q =行列式)或者通過它......? – HostileFork

+0

他們的方法對我來說沒有多大意義,因此很難適應我圖書館的運作方式。如果我理解了它背後的數學推理,我就可以很容易地適應它。我認爲這與我正在開始研究的部分旋轉有關。其他方法對我來說是有意義的,但這是我第一次無法理解它背後的方法。一般來說,只是好奇另一個大腦是否有'最好的方式'的想法。即使在發佈這個問題之後,我仍然非常關注它,並會在找到更好的方法時發佈我的代碼。 – user2927848

+0

[this](http://stackoverflow.com/questions/22479258/cholesky-decomposition-with-openmp)可能對你很有意思。 –

回答

3

算法通過高斯消去由n矩陣減少的n與上(或下)三角形式通常具有複雜度O(N^3) (其中^代表「權力」)。

還有其他的計算確定性的方法,例如評估特徵值集合(方陣矩陣的行列式等於它的特徵值的乘積)。對於一般的矩陣,完整的特徵值的計算也是 - 實際上 - O(n^3)。

然而,理論上,計算特徵值集的複雜度爲n^w,其中w在2和2.376之間 - 這意味着對於(很多)較大的矩陣,它將比使用高斯消元更快。看看James Demmel,Ioana Dumitriu和Olga Holtz在2007年11月第108期第1期第59-91頁的Numerische Mathematik的一篇文章「快速線性代數是否穩定」。如果Eigen使用複雜度較低的方法比O(n^3)更大的矩陣(我不知道,從來沒有理由調查這些事情),這將解釋你的觀察。

+0

謝謝這篇論文很有趣。我認爲'塊LU分解'是Eigen與8x8子塊一起使用的方法。那篇文章說時間複雜度大約是O(n^2.5)。我還將更多地考慮特徵值選項。 – user2927848

0

答案大多數地方似乎使用塊LU分解在相同的內存空間中創建下三角形和上三角形矩陣。根據您使用的塊的大小,它是〜O(n^2.5)。

下面是來自萊斯大學的一個功率點,它解釋了算法。

www.caam.rice.edu/~timwar/MA471F03/Lecture24.ppt

司由矩陣表示通過它的逆的乘法。

這個想法似乎是顯着地增加了n^2操作的數量,但是減少了數量m^3,這實際上降低了算法的複雜度,因爲m具有固定的小尺寸。

爲了有效地寫出這個問題,需要一點點寫,因爲要有效地完成這個任務需要「現場」算法,我還沒有寫過。

相關問題