固有値問題で MATLAB/Octave が C++ を使用してフロアを一掃するのはなぜですか?

タイトルの質問に対する答えが、私が愚かなことをしているということであることを願っています!

これが問題です。実対称行列のすべての固有値と固有ベクトルを計算したいと考えています。私は MATLAB (実際には、Octave を使用して実行しています) と C++ でコードを実装し、GNU Scientific Library を使用しています。両方の実装について、以下に完全なコードを提供します。

私が理解できる限り、GSL には BLAS API の独自の実装が付属しており (以下、これを GSLCBLAS と呼びます)、このライブラリを使用するために、以下を使用してコンパイルします:

g++ -O3 -lgsl -lgslcblas
 

GSL は、パフォーマンスを向上させるために、自己最適化 ATLAS ライブラリなどの代替 BLAS ライブラリを使用することをここで提案しています。 Ubuntu 12.04 を実行しており、Ubuntu リポジトリから ATLAS パッケージをインストールしました。この場合、以下を使用してコンパイルします:

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

3 つのケースすべてについて、ランダムに生成されたサイズ 100 から 1000 の行列を 100 刻みで使用して実験を行いました。サイズごとに、異なる行列を使用して 10 回の固有分解を実行し、所要時間を平均しました。結果は次のとおりです:

性能の違いはばかげています。サイズが 1000 の行列の場合、Octave は 1 秒未満で分解を実行します。 GSLCBLAS と ATLAS には約 25 秒かかります。

ATLAS ライブラリを間違って使用している疑いがあります。どんな説明でも大歓迎です。よろしくお願いします。

コードに関する注意事項:

• C++ 実装では、行列を対称にする必要はありません
  関数は行列の下三角部分のみを使用するためです
  。




• Octave では、行 triu(A) + triu(A, 1)' 行列が対称になるように強制します。




• 独自の Linux マシンで C++ コードをコンパイルする場合は、フラグ -lrt も追加する必要があります。 、 clock_gettime のため 関数。




• 残念ながら clock_gettime ではないと思います 他のプラットフォームでは終了します。 gettimeofday に変更することを検討してください .

オクターブ コード

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

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

いくつかのコードの回答

g++ -O3 -lgsl -lgslcblas 

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

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

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

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

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

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

(* 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]}]