Perché MATLAB/Octave pulisce il pavimento con C++ in Problemi autovalore?

Question

Spero che la risposta alla domanda nel titolo sia che sto facendo qualcosa di stupido!

Ecco il problema. Voglio calcolare tutti gli autovalori e gli autovettori di una matrice simmetrica reale. Ho implementato il codice in MATLAB (in realtà, lo eseguo usando Octave) e C++, usando lo GNU Scientific Library. Sto fornendo il mio codice completo qui sotto per entrambe le implementazioni.

Per quanto posso capire, GSL dotato di una propria implementazione delle API BLAS, (qui di seguito rimando a questo come GSLCBLAS) e di utilizzare questa libreria compilo usando:

g++ -O3 -lgsl -lgslcblas 

GSL suggerisce here utilizzare una libreria BLAS alternativa, ad esempio la libreria auto-ottimizzante ATLAS, per prestazioni migliorate. Sono in esecuzione Ubuntu 12.04 e ho installato i pacchetti ATLAS dal Ubuntu repository. In questo caso, compilo utilizzando:

esperimenti

g++ -O3 -lgsl -lcblas -latlas -lm 

Per tutti e tre i casi, sono eseguiti con matrici generate in modo casuale di dimensioni 100 a 1000 in passi di 100. Per ciascuna dimensione, effettuo 10 eigendecompositions con differenti matrici e media il tempo impiegato. I risultati sono questi:

La differenza di prestazioni è ridicolo. Per una matrice di dimensione 1000, Octave esegue la decomposizione in meno di un secondo; GSLCBLAS e ATLAS impiegano circa 25 secondi.

Ho il sospetto che possa utilizzare la libreria ATLAS in modo errato. Eventuali spiegazioni sono ben accette; Grazie in anticipo.

Alcune note sul codice:

• Nell'implementazione C++, non c'è bisogno di fare la matrice simmetrica , perché the function only uses the lower triangular part of it.

• In Octave, la riga triu(A) + triu(A, 1)' impone che la matrice sia simmetrica.

• Se si desidera compilare il codice C++ sulla propria macchina Linux, è inoltre necessario aggiungere la flag -lrt, a causa della funzione clock_gettime.

• Purtroppo non penso che le uscite clock_gettime su altre piattaforme. Valuta la possibilità di cambiarlo in gettimeofday.

codice ottava

K = 10; 

fileID = fopen('octave_out.txt','w'); 

for N = 100:100:1000 
    AverageTime = 0.0; 

    for k = 1:K 
     A = randn(N, N); 
     A = triu(A) + triu(A, 1)'; 
     tic; 
     eig(A); 
     AverageTime = AverageTime + toc/K; 
    end 

    disp([num2str(N), " ", num2str(AverageTime), "\n"]); 
    fprintf(fileID, '%d %f\n', N, AverageTime); 
end 

fclose(fileID); 

codice C++

#include <iostream> 
#include <fstream> 
#include <time.h> 

#include <gsl/gsl_rng.h> 
#include <gsl/gsl_randist.h> 
#include <gsl/gsl_eigen.h> 
#include <gsl/gsl_vector.h> 
#include <gsl/gsl_matrix.h> 

int main() 
{ 
    const int K = 10; 

    gsl_rng * RandomNumberGenerator = gsl_rng_alloc(gsl_rng_default); 
    gsl_rng_set(RandomNumberGenerator, 0); 

    std::ofstream OutputFile("atlas.txt", std::ios::trunc); 

    for (int N = 100; N <= 1000; N += 100) 
    { 
     gsl_matrix* A = gsl_matrix_alloc(N, N); 
     gsl_eigen_symmv_workspace* EigendecompositionWorkspace = gsl_eigen_symmv_alloc(N); 
     gsl_vector* Eigenvalues = gsl_vector_alloc(N); 
     gsl_matrix* Eigenvectors = gsl_matrix_alloc(N, N); 

     double AverageTime = 0.0; 
     for (int k = 0; k < K; k++) 
     { 
      for (int i = 0; i < N; i++) 
      { 
       for (int j = 0; j < N; j++) 
       { 
        gsl_matrix_set(A, i, j, gsl_ran_gaussian(RandomNumberGenerator, 1.0)); 
       } 
      } 

      timespec start, end; 
      clock_gettime(CLOCK_MONOTONIC_RAW, &start); 

      gsl_eigen_symmv(A, Eigenvalues, Eigenvectors, EigendecompositionWorkspace); 

      clock_gettime(CLOCK_MONOTONIC_RAW, &end); 
      double TimeElapsed = (double) ((1e9*end.tv_sec + end.tv_nsec) - (1e9*start.tv_sec + start.tv_nsec))/1.0e9; 
      AverageTime += TimeElapsed/K; 
      std::cout << "N = " << N << ", k = " << k << ", Time = " << TimeElapsed << std::endl; 
     } 
     OutputFile << N << " " << AverageTime << std::endl; 

     gsl_matrix_free(A); 
     gsl_eigen_symmv_free(EigendecompositionWorkspace); 
     gsl_vector_free(Eigenvalues); 
     gsl_matrix_free(Eigenvectors); 
    } 

    return 0; 
} 

Answer 1

Nell'implementazione C++, non v'è alcuna necessità di effettuare la matrice simmetrica , perché la funzione utilizza solo la parte triangolare inferiore di it.

Questo potrebbe non essere il caso. Nel reference, si precisa che:

int gsl_eigen_symmv (gsl_matrix * A, gsl_vector * eval, gsl_matrix * evec, gsl_eigen_symmv_workspace * w)

Questa funzione calcola gli autovalori e autovettori reale simmetrica matrice A. Spazio di lavoro aggiuntivo della dimensione appropriata deve essere fornito in w. La parte triangolare diagonale e inferiore di A viene distrutta durante il calcolo , ma la rigida parte triangolare superiore non viene referenziata. Gli autovalori sono memorizzati nella valutazione vettoriale e sono non ordinati. Gli autovettori corrispondenti sono memorizzati nelle colonne della matrice evec. Ad esempio, l'autovettore nella prima colonna corrisponde a il primo autovalore. Gli autovettori sono garantiti per essere reciprocamente ortogonali e normalizzati per unità di grandezza.

Sembra che sia necessario applicare un'operazione di simmetrizzazione simile in C++ per ottenere risultati almeno corretti sebbene sia possibile ottenere le stesse prestazioni.

Sul lato MATLAB, eigen valore decomposizione può essere più veloce grazie alla sua esecuzione multi-threaded come indicato in this reference:

incorporato Multithreading

algebra lineare e funzioni numeriche come fft, \ (mldivide), eig, svd e sort sono multithread in MATLAB. I calcoli con multithreading sono stati attivi per impostazione predefinita in MATLAB dalla release 2008a. Queste funzioni vengono eseguite automaticamente su più thread computazionali in una singola sessione MATLAB , consentendo loro di eseguire più rapidamente su macchine con tecnologia multicore . Inoltre, molte funzioni in Image Processing Toolbox ™ sono multithreaded.

Al fine di testare le prestazioni di MATLAB per single core, è possibile disattivare il multithreading per

File> Preferenze> Generale> Multithreading

in R2007a o più recente come indicato here.

Answer 2

Non sono d'accordo con il post precedente. Questo non è un problema di threading, questo è un problema di algoritmo. La ragione per cui matlab, R e ottava cancellano il pavimento con le librerie C++ è perché le loro librerie C++ usano algoritmi più complessi e migliori. Se leggete la pagina di ottava si può scoprire che cosa che fanno [1]:

Autovalori siano calcolati in diversi processo graduale che inizia con una decomposizione Hessenberg, seguiti da una decomposizione Schur, da cui gli autovalori sono evidenti. Gli autovettori, quando desiderato, sono calcolati da ulteriori manipolazioni della decomposizione di Schur.

La risoluzione dei problemi di autovalore/autovettore non è banale. In effetti, è una delle poche cose che "Ricette numeriche in C" ti consiglia di eseguire non. (P461). GSL è spesso lento, che è stata la mia prima risposta. ALGLIB è anche lento per la sua implementazione standard (sto ottenendo circa 12 secondi!):

#include <iostream> 
#include <iomanip> 
#include <ctime> 

#include <linalg.h> 

using std::cout; 
using std::setw; 
using std::endl; 

const int VERBOSE = false; 

int main(int argc, char** argv) 
{ 

    int size = 0; 
    if(argc != 2) { 
     cout << "Please provide a size of input" << endl; 
     return -1; 
    } else { 
     size = atoi(argv[1]); 
     cout << "Array Size: " << size << endl; 
    } 

    alglib::real_2d_array mat; 
    alglib::hqrndstate state; 
    alglib::hqrndrandomize(state); 
    mat.setlength(size, size); 
    for(int rr = 0 ; rr < mat.rows(); rr++) { 
     for(int cc = 0 ; cc < mat.cols(); cc++) { 
      mat[rr][cc] = mat[cc][rr] = alglib::hqrndnormal(state); 
     } 
    } 

    if(VERBOSE) { 
     cout << "Matrix: " << endl; 
     for(int rr = 0 ; rr < mat.rows(); rr++) { 
      for(int cc = 0 ; cc < mat.cols(); cc++) { 
       cout << setw(10) << mat[rr][cc]; 
      } 
      cout << endl; 
     } 
     cout << endl; 
    } 

    alglib::real_1d_array d; 
    alglib::real_2d_array z; 
    auto t = clock(); 
    alglib::smatrixevd(mat, mat.rows(), 1, 0, d, z); 
    t = clock() - t; 

    cout << (double)t/CLOCKS_PER_SEC << "s" << endl; 

    if(VERBOSE) { 
     for(int cc = 0 ; cc < mat.cols(); cc++) { 
      cout << "lambda: " << d[cc] << endl; 
      cout << "V: "; 
      for(int rr = 0 ; rr < mat.rows(); rr++) { 
       cout << setw(10) << z[rr][cc]; 
      } 
      cout << endl; 
     } 
    } 
} 

Se davvero bisogno di una libreria veloce, probabilmente bisogno di fare un po 'vera e propria caccia.

[1] http://www.gnu.org/software/octave/doc/interpreter/Basic-Matrix-Functions.html

Answer 3

Ho anche incontrato con il problema. La vera causa è che l'eig() in MATLAB non calcola gli autovettori, ma il codice di versione C sopra lo fa. Il diverso tempo trascorso può essere maggiore di un ordine di grandezza, come mostrato nella figura sottostante. Quindi il confronto non è giusto.

In Matlab, a seconda del valore restituito, la funzione effettiva chiamata sarà diversa. Per forzare il calcolo degli autovettori, è necessario utilizzare [V,D] = eig(A) (vedere i codici seguenti).

Il tempo effettivo per calcolare problema agli autovalori dipende fortemente dalle proprietà della matrice e dei risultati desiderati, come

• reale o complessa
• hermitiana/simmetrica o meno
• Dense o sparsi
• Solo autovalori, Autovettori, Autovalore massimo, ecc.
• Seriale o parallelo

Esistono algoritmi ottimizzati per ciascuno dei casi precedenti. Nel file gsl, questi algoritmi sono picked manually, pertanto una selezione errata ridurrà significativamente le prestazioni. Qualche classe di wrapper C++ o un linguaggio come matlab e mathematica sceglieranno la versione ottimizzata attraverso alcuni metodi.

Inoltre, Matlab e Mathematica hanno utilizzato la parallelizzazione. Questi sono ulteriormente ampliare il divario che vedi poche volte, a seconda della macchina. È ragionevole dire che il calcolo di autovalori e autovettori di un complesso generale 1000x1000 è di circa un secondo e dieci secondi, senza parallelizzazione.

Fig. Confronta Matlab e C. "+ vec" indica che i codici includevano i calcoli degli autovettori. La CPU% è l'osservazione molto approssimativa dell'uso della CPU a N = 1000, che è superiore del 800%, sebbene si supponga che utilizzino completamente tutti gli 8 core. Il divario tra Matlab e C è inferiore a 8 volte.

Fig. Confrontare il tipo di matrice diverso in Mathematica. Algoritmi selezionati automaticamente dal programma.

Matlab (con il calcolo di autovettori)

K = 10; 

fileID = fopen('octave_out.txt','w'); 

for N = 100:100:1000 
    AverageTime = 0.0; 

    for k = 1:K 
     A = randn(N, N); 
     A = triu(A) + triu(A, 1)'; 
     tic; 
     [V,D] = eig(A); 
     AverageTime = AverageTime + toc/K; 
    end 

    disp([num2str(N), ' ', num2str(AverageTime), '\n']); 
    fprintf(fileID, '%d %f\n', N, AverageTime); 
end 

fclose(fileID); 

C++ (senza il calcolo degli autovettori)

#include <iostream> 
#include <fstream> 
#include <time.h> 

#include <gsl/gsl_rng.h> 
#include <gsl/gsl_randist.h> 
#include <gsl/gsl_eigen.h> 
#include <gsl/gsl_vector.h> 
#include <gsl/gsl_matrix.h> 

int main() 
{ 
    const int K = 10; 

    gsl_rng * RandomNumberGenerator = gsl_rng_alloc(gsl_rng_default); 
    gsl_rng_set(RandomNumberGenerator, 0); 

    std::ofstream OutputFile("atlas.txt", std::ios::trunc); 

    for (int N = 100; N <= 1000; N += 100) 
    { 
     gsl_matrix* A = gsl_matrix_alloc(N, N); 
     gsl_eigen_symm_workspace* EigendecompositionWorkspace = gsl_eigen_symm_alloc(N); 
     gsl_vector* Eigenvalues = gsl_vector_alloc(N); 

     double AverageTime = 0.0; 
     for (int k = 0; k < K; k++) 
     { 
      for (int i = 0; i < N; i++) 
      { 
       for (int j = i; j < N; j++) 
       { 
        double rn = gsl_ran_gaussian(RandomNumberGenerator, 1.0); 
        gsl_matrix_set(A, i, j, rn); 
        gsl_matrix_set(A, j, i, rn); 
       } 
      } 

      timespec start, end; 
      clock_gettime(CLOCK_MONOTONIC_RAW, &start); 

      gsl_eigen_symm(A, Eigenvalues, EigendecompositionWorkspace); 

      clock_gettime(CLOCK_MONOTONIC_RAW, &end); 
      double TimeElapsed = (double) ((1e9*end.tv_sec + end.tv_nsec) - (1e9*start.tv_sec + start.tv_nsec))/1.0e9; 
      AverageTime += TimeElapsed/K; 
      std::cout << "N = " << N << ", k = " << k << ", Time = " << TimeElapsed << std::endl; 
     } 
     OutputFile << N << " " << AverageTime << std::endl; 

     gsl_matrix_free(A); 
     gsl_eigen_symm_free(EigendecompositionWorkspace); 
     gsl_vector_free(Eigenvalues); 
    } 

    return 0; 
} 

Mathematica

(* Symmetric real matrix + eigenvectors *) 
Table[{NN, Mean[Table[(
    M = Table[Random[], {i, NN}, {j, NN}]; 
    M = M + Transpose[Conjugate[M]]; 
    AbsoluteTiming[Eigensystem[M]][[1]] 
    ), {K, 10}]] 
    }, {NN, Range[100, 1000, 100]}] 

(* Symmetric real matrix *) 
Table[{NN, Mean[Table[(
    M = Table[Random[], {i, NN}, {j, NN}]; 
    M = M + Transpose[Conjugate[M]]; 
    AbsoluteTiming[Eigenvalues[M]][[1]] 
    ), {K, 10}]] 
    }, {NN, Range[100, 1000, 100]}] 

(* Asymmetric real matrix *) 
Table[{NN, Mean[Table[(
    M = Table[Random[], {i, NN}, {j, NN}]; 
    AbsoluteTiming[Eigenvalues[M]][[1]] 
    ), {K, 10}]] 
    }, {NN, Range[100, 1000, 100]}] 

(* Hermitian matrix *) 
Table[{NN, Mean[Table[(
    M = Table[Random[] + I Random[], {i, NN}, {j, NN}]; 
    M = M + Transpose[Conjugate[M]]; 
    AbsoluteTiming[Eigenvalues[M]][[1]] 
    ), {K, 10}]] 
    }, {NN, Range[100, 1000, 100]}] 

(* Random complex matrix *) 
Table[{NN, Mean[Table[(
    M = Table[Random[] + I Random[], {i, NN}, {j, NN}]; 
    AbsoluteTiming[Eigenvalues[M]][[1]] 
    ), {K, 10}]] 
    }, {NN, Range[100, 1000, 100]}]