2016-03-28 41 views
3

Ich schreibe eine Bibliothek, wo ich einige grundlegende NxN-Matrix-Funktionalität haben möchte, die keine Abhängigkeiten hat und es ist ein bisschen ein Lernprojekt. Ich vergleiche meine Leistung mit Eigen. Ich konnte in der Lage sein, ziemlich gleich zu sein und sogar seine Leistung auf einer Paarfront mit SSE2 zu übertreffen und mit AVX2 an einigen Fronten zu schlagen (es benutzt nur SSE2, also nicht super überraschend).Schnellster Weg Determinante zu berechnen?

Mein Problem ist, ich Gaußsche Eliminations bin mit einem Ober diagonalisierte Matrix dann die Diagonale Multiplikation erstellen die determinant.I schlagen Eigen für N < 300 zu bekommen, aber danach bläst Eigen mich weg und es wird nur noch schlimmer als die Matrizen größer werden. Angesichts der Tatsache, dass auf den gesamten Speicher sequenziell zugegriffen wird und die Zerlegung des Compilers nicht schrecklich aussieht, glaube ich nicht, dass es sich um ein Optimierungsproblem handelt. Es gibt mehr Optimierung, die durchgeführt werden kann, aber die Timings sehen viel mehr wie ein algorithmisches Timing-Komplexitäts-Problem aus oder es gibt einen großen SSE-Vorteil, den ich nicht sehe. Einfach die Loops ein wenig abzurollen hat mir nicht viel gebracht.

Gibt es einen besseren Algorithmus zur Berechnung von Determinanten?

Scalar Code

/* 
    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

Vorausgesetzt, Sie sind klar bereit, in Details zu graben, und Eigen ist Open Source ... [Warum nicht schauen, was es tut] (https://github.com/RLovelett/eigen/search?utf8 =% E2% 9C% 93 & q = Determinante) oder durchschreiten ...? – HostileFork

+0

Ihre Methode macht für mich keinen Sinn und macht es schwierig, mich an die Arbeitsweise meiner Bibliothek anzupassen. Wenn ich die mathematische Argumentation dahinter verstehen würde, wäre ich in der Lage, sie einfach genug anzupassen. Ich denke, es hat etwas mit teilweisem Schwenken zu tun, in das ich zu schauen beginne. Die anderen Methoden haben für mich Sinn ergeben, aber dies ist das erste Mal, dass ich die Methodik dahinter nicht verstehen konnte. Im Allgemeinen, nur neugierig, ob ein anderes Gehirn eine Vorstellung von "dem besten Weg" hat. Auch nachdem ich diese Frage gepostet habe, schaue ich immer noch sehr genau hin und poste meinen Code, wenn ich einen besseren Weg finde. – user2927848

+0

[dies] (http://stackoverflow.com/questions/22479258/cholesky-decomposition-with-openmp) könnte für Sie interessant sein. –

Antwort

3

Algorithms ein n durch n Matrix obere (oder untere) Dreiecksform durch Gaußsche Eliminations allgemeine Komplexität von O hat zu reduzieren (n^3) (wobei^für "Macht" steht).

Es gibt alternative Ansätze zum Berechnen von bestimmten, z. B. das Bewerten der Menge von Eigenwerten (die Determinante einer quadratischen Matrix ist gleich dem Produkt ihrer Eigenwerte). Für allgemeine Matrizen ist die Berechnung des vollständigen Satzes von Eigenwerten auch - praktisch - O (n^3).

In der Theorie hat die Berechnung der Menge der Eigenwerte jedoch eine Komplexität von, wobei w zwischen 2 und 2.376 liegt - was bedeutet, dass sie für (viel) größere Matrizen schneller ist als die Gaußsche Eliminierung. Siehe dazu den Artikel "Schnelle lineare Algebra ist stabil" von James Demmel, Ioana Dumitriu und Olga Holtz in Numerische Mathematik, Band 108, Heft 1, S. 59-91, November 2007. Wenn Eigen einen Ansatz mit weniger Komplexität verwendet als O (n^3) für größere Matrizen (ich weiß nicht, nie Grund gehabt, solche Dinge zu untersuchen), die Ihre Beobachtungen erklären würden.

+0

Danke, dass Papier interessant ist. Ich denke, Block LU Factorization ist die Methode, die Eigen mit 8x8 Subblöcken verwendet. In diesem Artikel steht, dass die Komplexität der Zeit bei O (n^2,5) liegt. Ich werde auch mehr auf die Eigenwertoption eingehen. – user2927848

0

Die Antwort die meisten Orte scheinen Block LU Factorization zu verwenden, um eine untere Dreiecksmatrix und eine obere Dreiecksmatrix im selben Speicherplatz zu erstellen. Es ist ~ O (n^2,5) abhängig von der Größe des verwendeten Blocks.

Hier ist ein Power Point von Rice University, der den Algorithmus erklärt.

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

Division durch eine Matrixeinrichtung durch ihre inverse Multiplikation.

Die Idee scheint zu sein, die Anzahl von n^2 Operationen signifikant zu erhöhen, aber die Anzahl m^3 zu verringern, was in der Tat die Komplexität des Algorithmus verringert, da m eine feste kleine Größe hat.

Gehen Sie ein wenig zu nehmen, um dies auf eine effiziente Art und Weise zu schreiben, da es effizient macht 'in Place' Algorithmen, die ich noch nicht geschrieben habe.

Verwandte Themen