Compute matrix eigenvalues on gsl

Question

I have define a mateig to compute matrix eigenvalues on 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;
}

I have test two matrix with the output (aka ref_eigQ[]) of matlab function eig(). for example

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

But it did not work with this matrix

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

But mateig() return [ 5.701562, -0.701562,-2]. I do not know why. Plz help me. Thanks.

Answer 1

Well, I get it.

The problem is the function gsl_eigen_symmv(), which is only for symmetric matrices.

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

While, matrix Q is not a symmetric matrix. But the function gsl_eigen_symmv() will process it as a symmetric matrix. It automatically convert martrix Q into a symmetric form:

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

It lead the error. gsl_eigen_nonsymmv() will be ok for matrix Q.