在 gsl 上计算矩阵特征值

Question

我定义了一个mateig来计算gsl

上的矩阵特征值

#include <gsl/gsl_math.h>                                                                                                                                                                                        #include <gsl/gsl_sf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
/*
 * @brief matlab eig()
 *
 */
static int mateig(const double *Q, int n, int m, double *eigQ){
    int i,j;
    gsl_matrix *Qm = gsl_matrix_calloc(n, m);
    gsl_vector_view eval = gsl_vector_view_array (eigQ,n);
    gsl_matrix *evec = gsl_matrix_alloc (n,m);

    for(i=0;i<n;i++){
        for(j=0;j<m;j++){
            gsl_matrix_set(Qm, i, j,  Q[i*n + j]/1.0);
        }
    }

    gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (n);
    gsl_eigen_symmv (Qm, &eval.vector, evec, w);
    gsl_eigen_symmv_free (w);


    gsl_eigen_symmv_sort (&eval.vector, evec, GSL_EIGEN_SORT_VAL_DESC);

    gsl_matrix_free (Qm);
    gsl_matrix_free (evec);
    return 1;
}

我用 matlab 函数 eig() 的输出（又名 ref_eigQ[]）测试了两个矩阵。例如

    double Q[] = {
        1.0 , 1/2.0, 1/3.0, 1/4.0,
        1/2.0, 1/3.0, 1/4.0, 1/5.0,
        1/3.0, 1/4.0, 1/5.0, 1/6.0,
        1/4.0, 1/5.0, 1/6.0, 1/7.0 };
    double ref_eigQ[]={ 1.500214, 0.169141, 0.006738, 0.000097 };

    double Q[] = {
        1.0000,    0.5000 ,   0.3333  ,  0.2500,
        0.5000,    1.0000 ,   0.6667  ,  0.5000,
        0.3333,    0.6667 ,   1.0000  ,  0.7500,
        0.2500,    0.5000 ,   0.7500  ,  1.0000
    };
    double ref_eigQ[]={
        2.536162474486201,
        0.848229155477913,
        0.407832884117875,
        0.207775485918012
    };

但它不适用于此矩阵

    double Q[] = {
        1,2,3,
        2,1,3,
        3,2,1}; //det = 12
    double ref_eigQ[]={ 6,-1,-2 };

但是mateig()return[ 5.701562, -0.701562,-2]。我不知道为什么。请帮助我。谢谢。

Answer 1

嗯，我明白了。

问题是函数 gsl_eigen_symmv()，它仅适用于对称矩阵。

 double Q[] = {
        1,2,3,
        2,1,3,
        3,2,1};

而矩阵Q不是对称矩阵。但是函数 gsl_eigen_symmv() 会将其处理为对称矩阵。它自动将矩阵 Q 转换为对称形式：

 double Q[] = {
        1,2,3,
        2,1,2,
        3,2,1};

导致错误。 gsl_eigen_nonsymmv() 可以用于矩阵 Q。

在 gsl 上计算矩阵特征值

Compute matrix eigenvalues on gsl

c

gsl