我希望标题中问题的答案是我在做一些愚蠢的事情!
这就是问题所在。我想计算一个实数对称矩阵的所有特征值和特征向量。我已经在 MATLAB(实际上,我使用 Octave 运行它)和 C++ 中实现了代码,使用 GNU Scientific Library .我在下面为这两种实现提供了我的完整代码。
据我所知,GSL 带有自己的 BLAS API 实现(以下我将其称为 GSLCBLAS)并使用我编译的这个库:
g++ -O3 -lgsl -lgslcblas
GSL 建议 here使用替代的 BLAS 库,例如自优化 ATLAS库,以提高性能。我正在运行 Ubuntu 12.04,并且已经从 Ubuntu repository 安装了 ATLAS 包。 .在这种情况下,我编译使用:
g++ -O3 -lgsl -lcblas -latlas -lm
对于所有这三种情况,我都以 100 的步长对大小为 100 到 1000 的随机生成的矩阵进行了实验。对于每个大小,我使用不同的矩阵执行 10 次特征分解,并平均花费的时间。结果如下:
性能上的差异是荒谬的。对于大小为 1000 的矩阵,Octave 在不到一秒的时间内执行分解; GSLCBLAS 和 ATLAS 大约需要 25 秒。
我怀疑我可能错误地使用了 ATLAS 库。欢迎任何解释;提前致谢。
代码的一些注释:
在C++实现中,不需要制作矩阵 对称,因为 the function only uses the lower triangular part of it .
在 Octave 中,
triu(A) + triu(A, 1)'
行强制矩阵是对称的。如果您想在自己的 Linux 机器上编译 C++ 代码,您还需要添加标志
-lrt
,因为有clock_gettime
功能。不幸的是,我认为
clock_gettime
不会在其他平台上退出。考虑将其更改为gettimeofday
。
Octave 码
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;
}
最佳答案
我不同意之前的帖子。这不是线程问题,这是算法问题。 matlab、R 和 octave 使用 C++ 库的原因是因为它们的 C++ 库使用更复杂、更好的算法。如果您阅读 octave 页面,您会发现 他们 做了什么[1]:
Eigenvalues are computed in a several step process which begins with a Hessenberg decomposition, followed by a Schur decomposition, from which the eigenvalues are apparent. The eigenvectors, when desired, are computed by further manipulations of the Schur decomposition.
解决特征值/特征向量问题并非易事。事实上,它是“C 中的数值食谱”建议您不要自己实现的少数几件事之一。 (第 461 页)。 GSL 通常很慢,这是我最初的 react 。 ALGLIB 的标准实现也很慢(我大约需要 12 秒!):
#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;
}
}
}
如果你真的需要一个快速的图书馆,可能需要做一些真正的狩猎。
[1] http://www.gnu.org/software/octave/doc/interpreter/Basic-Matrix-Functions.html
关于c++ - 为什么 MATLAB/Octave 在特征值问题中用 C++ 擦地板?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/18009056/