GSL特征值顺序

Question

我正在使用GSL库中的函数gsl_eigen_nonsymm和/或gsl_eigen_symm来找到L x L矩阵的特征值,该矩阵M[i][j]也是时间的函数,t = 1,....,N因此我必须M[i][j][t]得到每个ti的特征值,分配一个L x L矩阵E[i][j] = M[i][j][t]并对其进行对角化对于每一个.

问题是程序在一些迭代后给出了不同顺序的特征值.例如(L = 3),如果在t = 0我得到eigen[t = 0] = {l1,l2,l3}(0)在 t = 1i可以得到eigen[t = 1] = {l3,l2,l1}(1)而我需要总是有{l1,l2,l3}(t) 更具体的:考虑矩阵M (t) ) = {{0,t,t},{t,0,2t},{t,t,0}}的特征值将始终是(approximatevly)l1 = -1.3 t , l2 = -t , l3 = 2.3 t当我试图对角化它(与下面的代码)我在特征值的结果中多次交换.有没有办法防止它？我不能只按它的大小对它们进行排序我需要它们始终处于相同的顺序(无论它是什么)先验.(下面的代码只是解释我的问题的一个例子)

编辑:我不能只对它们进行排序,因为先验我不知道它们的价值,也不是l1<l2<l3因为统计波动它们是否可靠地具有类似的结构,这就是为什么我想知道是否有办法制作算法行为总是以相同的方式运行,以便特征值的顺序始终相同,或者是否有一些技巧可以使它发生.

为了更清楚,我将尝试重新描述我在这里提出的玩具问题.我们有一个矩阵,它取决于时间,我可能天真地,期望得到 lambda_1(t).....lambda_N(t),而我所看到的是算法经常在不同的时间交换特征值,所以如果t = 1 I've got ( lambda_1,lambda_2,lambda_3 )(1) at time t = 2 (lambda_2,lambda_1,lambda_3)(2)这样,例如我想看看lambda_1如何进化时间我不能,因为算法在不同的时间混合特征值.下面的程序只是我问题的分析玩具示例:下面矩阵的特征值是l1 = -1.3 t , l2 = -t , l3 = 2.3 t但是程序可能会给我输出(-1.3,-1,2.3)(1), (-2,-2.6,4.6)(2), etc如前所述,我想知道是否有办法使程序命令特征值尽管它们的实际数值总是以相同的方式,所以我总是得到(l1,l2,l3)组合.我希望现在更清楚,请告诉我是不是.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sort_vector.h>

main() {
    int L = 3, i, j, t;
    int N = 10;
    double M[L][L][N];
    gsl_matrix *E = gsl_matrix_alloc(L, L);
    gsl_vector_complex *eigen = gsl_vector_complex_alloc(L);
    gsl_eigen_nonsymm_workspace *  w = gsl_eigen_nonsymm_alloc(L);

    for(t = 1; t <= N; t++) {
        M[0][0][t-1] = 0;
        M[0][1][t-1] = t;
        M[0][2][t-1] = t;
        M[1][0][t-1] = t;
        M[1][1][t-1] = 0;
        M[1][2][t-1] = 2.0 * t;
        M[2][1][t-1] = t;
        M[2][0][t-1] = t;
        M[2][2][t-1] = 0;

        for(i = 0; i < L; i++) {
            for(j = 0; j < L; j++) {
                gsl_matrix_set(E, i, j, M[i][j][t - 1]);
            }
        }

        gsl_eigen_nonsymm(E, eigen, w); /*diagonalize E which is M at t fixed*/
        printf("#%d\n\n", t);

        for(i = 0; i < L; i++) {
            printf("%d\t%lf\n", i, GSL_REAL(gsl_vector_complex_get(eigen, i)))
        }
        printf("\n");
   }
}

Answer 1

我不认为你能做你想做的事。随着t变化，输出也会变化。

我原来的答案提到了指针的排序，但查看数据结构并没有帮助。计算出特征值后，这些值将存储在 E 中。您可以如下查看它们。

gsl_eigen_nonsymm(E, eigen, w);
double *mdata = (double*)E->data;
printf("mdata[%i]    \t%lf\n", 0, mdata[0]);
printf("mdata[%i]    \t%lf\n", 4, mdata[4]);
printf("mdata[%i]    \t%lf\n", 8, mdata[8]);

以下代码是特征向量中的数据的布局方式......

double *data  = (double*)eigen->data;
for(i = 0; i < L; i++) {
  printf("%d n  \t%zu\n", i, eigen->size);
  printf("%d    \t%lf\n", i, GSL_REAL(gsl_vector_complex_get(eigen, i)));
  printf("%d r  \t%lf\n", i, data[0]);
  printf("%d i  \t%lf\n", i, data[1]);
  printf("%d r  \t%lf\n", i, data[2]);
  printf("%d i  \t%lf\n", i, data[3]);
  printf("%d r  \t%lf\n", i, data[4]);
  printf("%d i  \t%lf\n", i, data[5]);
}   

如果，并且当您看到顺序更改、mdata更改中的数据顺序以及更改中的顺序时，您可以检查这一点data，那么该算法没有固定的顺序，即您无法执行您要求它执行的操作。如果顺序没有改变mdata并且发生了变化，data那么你就有了一个解决方案，但我真的怀疑情况会是这样。