为什么MATLAB / Octave在特征值问题中用C ++擦拭地板？
我希望题目中问题的答案是我正在做一些愚蠢的事情！
这是问题。 我想计算一个真正的对称矩阵的所有特征值和特征向量。 我已经在MATLAB中实现了代码（实际上，我使用Octave来运行它）以及使用GNU科学库的C ++。 我在下面提供了两个实现的完整代码。
据我所知，GSL自带BLAS API实现（以下我称之为GSLCBLAS），并使用我编译的库：
g++ -O3 -lgsl -lgslcblas
GSL建议在此使用替代的BLAS库，例如自我优化的ATLAS库，以提高性能。 我正在运行Ubuntu 12.04，并且已经从Ubuntu存储库安装了ATLAS软件包。 在这种情况下，我编译使用：
g++ -O3 -lgsl -lcblas -latlas -lm
对于所有三种情况，我都用100到100的随机生成矩阵进行了实验，步长为100。对于每个大小，我用不同的矩阵执行10次特征分解，并平均所花费的时间。 结果如下：
表现的差异是荒谬的。 对于大小为1000的矩阵，Octave在一秒之内执行分解; GSLCBLAS和ATLAS大约需要25秒。
我怀疑我可能正在错误地使用ATLAS库。 任何解释都是值得欢迎的。 提前致谢。
关于代码的一些说明：
在C ++实现中，不需要使矩阵对称，因为函数只使用它的下三角部分。
在八度中，行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 %fn', 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 ++库使用更复杂，更好的算法。 如果你阅读八度页面，你可以找出他们做什么[1]：
特征值是在几步过程中计算出来的，该过程从Hessenberg分解开始，然后是Schur分解，从中可以看出特征值。 当需要时，特征向量通过对Schur分解的进一步处理来计算。
解决特征值/特征向量问题并不重要。 实际上，它是“C中的数字食谱”中少数几件事之一，建议您不要实现自己。 （P461）。 GSL通常很慢，这是我最初的回应。 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
我也遇到过这个问题。 真正的原因是matlab中的eig（）不会计算特征向量，但是上面的C版本代码不会。 如下图所示，花费的时间不同可能大于一个数量级。 所以比较是不公平的。
在Matlab中，根据返回值，调用的实际函数将会不同。 为了强制计算特征向量，应该使用[V，D] = eig（A）（见下面的代码）。
计算特征值问题的实际时间在很大程度上取决于矩阵特性和期望的结果，例如
有针对上述每种情况优化的算法。 在gsl中，这些算法是手动挑选的，所以错误的选择会显着降低性能。 一些C ++包装类或诸如matlab和mathematica之类的语言将通过一些方法来选择优化版本。
另外，Matlab和Mathematica也使用了并行。 根据机器的不同，这会进一步扩大您看到的几次间隔。 可以合理地说，一般复数1000x1000的特征值和特征向量的计算大约是一秒和十秒，没有并行化。
图。比较Matlab和C.“+ vec”表示代码包含特征向量的计算。 CPU％是N = 1000时CPU使用率的非常粗略的观察结果，其上限为800％，尽管它们应该完全使用全部8个内核。 Matlab和C之间的差距小于8倍。
图。比较Mathematica中不同的矩阵类型。 算法由程序自动选择。
Matlab（用特征向量的计算）
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 %fn', 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_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]}]
在C ++实现中，不需要使矩阵对称，因为函数只使用它的下三角部分。
情况可能并非如此。 在参考文献中指出：
int gsl_eigen_symmv（gsl_matrix * A，gsl_vector * eval，gsl_matrix * evec，gsl_eigen_symmv_workspace * w）
该函数计算实对称矩阵A的特征值和特征向量。 必须在w中提供适当大小的附加工作空间。 A的对角线和下三角部分在计算过程中被破坏，但严格的上三角部分未被引用。 特征值存储在矢量eval中，并且是无序的。 相应的特征向量存储在矩阵平面的列中。 例如，第一列中的特征向量对应于第一特征值。 特征向量保证相互正交并归一化为单位量级。
看起来您也需要在C ++中应用类似的对称化操作，以获得至少正确的结果，尽管您可以获得相同的性能。
在MATLAB方面，由于其多线程执行特征值分解可能会更快，如本参考中所述：
内置多线程
线性代数和数值函数，如fft，（mldivide），eig，svd和sort在MATLAB中是多线程的。 自2008年发布以来，MATLAB中默认使用多线程计算。 这些函数可以在单个MATLAB会话中的多个计算线程上自动执行，从而可以在启用多核的机器上更快地执行。 此外，Image Processing Toolbox™中的许多功能都是多线程的。
为了测试MATLAB单核的性能，可以通过禁用多线程
文件>首选项>常规>多线程
在R2007a或更新如此处所述。
链接地址: http://www.djcxy.com/p/85967.html上一篇: Why does MATLAB/Octave wipe the floor with C++ in Eigenvalue Problems?