为什么MATLAB/Octave在特征值问题中用C++擦拭地板？

Question

我希望标题中的问题的答案是我做了一些愚蠢的事情!

这是问题所在.我想计算一个真实对称矩阵的所有特征值和特征向量.我已经使用GNU Scientific Library在MATLAB中实现了代码(实际上,我使用Octave运行它)和C++ .我在下面提供了完整的代码用于两种实现.

据我所知,GSL带有自己的BLAS API实现,(以下我将其称为GSLCBLAS)并使用我使用以下编译的库:

g++ -O3 -lgsl -lgslcblas

GSL表明这里使用替代BLAS库,如自我优化ATLAS库,以提高性能.我正在运行Ubuntu 12.04,并已从Ubuntu存储库安装了ATLAS软件包.在这种情况下,我编译使用:

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

对于所有三种情况,我已经使用随机生成的大小为100到1000的矩阵进行了实验,步长为100.对于每个大小,我执行10个具有不同矩阵的特征分解,并平均所花费的时间.结果如下:

性能上的差异是荒谬的.对于大小为1000的矩阵,Octave在一秒钟内执行分解; GSLCBLAS和ATLAS大约需要25秒.

我怀疑我可能错误地使用了ATLAS库.欢迎任何解释; 提前致谢.

关于代码的一些注意事项:

• 在C++实现中,不需要使矩阵对称,因为该函数仅使用它的下三角部分.

• 在Octave中,该行triu(A) + triu(A, 1)'强制矩阵是对称的.

• 如果您希望将C++代码编译为您自己的Linux机器,则还需要添加该标志-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;
}

Answer 1

我不同意上一篇文章.这不是一个线程问题,这是一个算法问题.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

Answer 2

\n

在 C++ 实现中，不需要使矩阵对称，因为该函数仅使用 的下三角部分。

\n

\n

情况可能并非如此。在参考文献中，指出：

\n

\n

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

\n

此函数计算实对称矩阵\nA的特征值和特征向量。必须在 w 中提供适当大小的附加工作空间。\nA 的对角线和下三角部分在计算过程中被破坏，\n但严格的上三角部分不被引用。\n特征值存储在向量 eval 中并且是无序的。相应的特征向量存储在矩阵的列中。例如，第一列中的特征向量对应于第一个特征值。保证特征向量相互正交并归一化为单位幅度。

\n

\n

看来您还需要在 C++ 中应用类似的对称化操作，以便至少获得正确的结果，尽管您可以获得相同的性能。

\n

在 MATLAB 方面，特征值分解可能会更快，因为它是多线程执行的，如本参考文献中所述：

\n

\n

内置多线程

\n

线性代数和数值函数（例如 fft、\\ (mldivide)、eig、\nsvd 和 sort）在 MATLAB 中是多线程的。自 2008a 版本以来，MATLAB\n默认启用多线程计算。这些函数在单个 MATLAB 会话中的多个计算线程上自动执行，从而使它们能够在支持多核的计算机上更快地执行。此外，Image\nProcessing Toolbox\xe2\x84\xa2 中的许多函数都是多线程的。

\n

\n

为了测试 MATLAB 单核的性能，您可以通过以下方式禁用多线程：

\n

\n

文件>首选项>常规>多线程

\n

\n

在 R2007a 或更高版本中，如此处所述。

\n