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] 并为每个 t 对角化它。

问题是程序在一些迭代之后以不同的顺序给出特征值。例如（L = 3）如果在t = 0 我得到eigen[t = 0] = {l1,l2,l3}(0) 在t = 1 我可能会得到eigen[t = 1] = {l3,l2,l1}(1) 而我需要总是有{l1,l2,l3}(t) 更具体地说：考虑矩阵M (t) ) = {{0,t,t},{t,0,2t},{t,t,0}}，特征值将始终（近似）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

根据文档，这些函数返回无序：

https://www.gnu.org/software/gsl/manual/html_node/Real-Symmetric-Matrices.html

此函数计算实对称矩阵 A 的特征值。必须在 w 中提供适当大小的附加工作空间。 A的对角线和下三角部分在计算过程中被破坏，但严格的上三角部分没有被引用。特征值存储在向量 eval 中，并且是无序的。

即使是返回有序结果的函数，也可以通过简单的升序/降序来实现：

https://www.gnu.org/software/gsl/manual/html_node/Sorting-Eigenvalues-and-Eigenvectors.html

该函数同时将向量eval中存储的特征值和矩阵evec各列中存储的对应实特征向量根据参数sort_type的值进行升序或降序排序，如上图所示。

如果您正在寻找特征值的时间演化，只需像您一直在做的那样解决与时间相关的表示，例如：

lambda_1(t).....lambda_N(t)

对于您的简单时间标量示例，

l1 = -1.3 t , l2 = -t , l3 = 2.3 t

实际上，您已经对所有可能的解决方案进行了参数化，并且因为您已经为它们分配了标识符ln，所以您不会遇到退化问题。即使任何M[i][j] 是 t 的非线性函数，也没有关系，因为系统本身是线性的，并且解决方案纯粹由特征方程计算（在求解 lambda 时将保持 t 不变）。

Answer 2

您的问题零意义。特征值对它们没有任何内在顺序。在我看来，您想定义 M_t 的特征值类似于 L_1(M_t),..., L_n(M_t)，然后跟踪它们如何随时间变化。假设您驱动 M_t 的过程是连续的，那么您的特征值也将是连续的。换句话说，当您对 M_t 进行小的更改时，它们不会发生显着变化。因此，如果您通过强制 L_1

这是跟踪特征向量的另一种方法，可能会更好。为此，假设您的特征向量是 v_i，具有组件 v_ij 。你要做的是首先“标准化”你的特征向量，使得 v_i1 是非负的，即适当地翻转每个特征向量的符号。这将通过对每个特征向量的第一个分量 v_i1 的排序来定义特征值的排序。这样，您仍然可以跟踪相互交叉的特征值。但是，如果您的特征向量在第一个组件上交叉，您就有麻烦了。

Answer 3

我不认为你可以做你想做的事。随着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 中的订单发生了变化，那么您有一个解决方案，但我真的怀疑会是这种情况。