GSL 特征值顺序
GSL eigenvalues order
我正在使用 GSL 库中的函数 gsl_eigen_nonsymm and/or gsl_eigen_symm 来查找 L x L 矩阵 M[i][j]
的特征值,它也是时间 t = 1,....,N
所以我有 M[i][j][t]
来获取每个 t 的特征值 我分配一个 L x L 矩阵 E[i][j] = M[i][j][t]
并为每个 t 对角化它。
问题是程序在一些迭代后以不同的顺序给出了特征值。例如 (L = 3) 如果在 t = 0
我在 t = 1
得到 eigen[t = 0] = {l1,l2,l3}(0)
我可能会得到 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");
}
}
我不认为你可以为所欲为。随着 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
中改变,那么你有一个解决方案,但我真的怀疑情况是否如此。
你的问题毫无意义。特征值对它们没有任何固有的顺序。在我看来,您想定义 M_t 的特征值类似于 L_1(M_t),..., L_n(M_t) 然后跟踪它们如何随时间变化。假设您的过程驱动 M_t 是连续的,那么您的特征值也是连续的。换句话说,当您对 M_t 进行小的更改时,它们不会发生显着变化。因此,如果您通过强制执行 L_1 < L_2... < L_n 来定义排序,那么此排序不会因 t 的微小变化而改变。当您有两个特征值交叉时,您需要决定如何分配更改。如果 "random fluctuations" 大于特征值之间的典型距离,那么这基本上是不可能的。
这是另一种跟踪特征向量的方法,可能会更好。为此,假设您的特征向量为 v_i,分量为 v_ij 。你首先要做的是 "normalize" 你的特征向量使得 v_i1 是非负的,即适当地翻转每个特征向量的符号。这将通过对 v_i1(每个特征向量的第一个分量)的排序来定义特征值的排序。这样您仍然可以跟踪相互交叉的特征值。但是,如果您的特征向量在第一个组件上交叉,您就有麻烦了。
根据文档,这些函数 return 无序:
https://www.gnu.org/software/gsl/manual/html_node/Real-Symmetric-Matrices.html
This function computes the eigenvalues of the real symmetric matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The eigenvalues are stored in the vector eval and are unordered.
即使是 return 排序的函数,也可以通过简单的 ascending/descending 量级来实现:
https://www.gnu.org/software/gsl/manual/html_node/Sorting-Eigenvalues-and-Eigenvectors.html
This function simultaneously sorts the eigenvalues stored in the vector eval and the corresponding real eigenvectors stored in the columns of the matrix evec into ascending or descending order according to the value of the parameter sort_type as shown above.
如果您正在寻找特征值的时间演化,就像您一直在做的那样并求解时间相关的表示,例如:
lambda_1(t).....lambda_N(t)
对于简单的时间标量示例,
l1 = -1.3 t , l2 = -t , l3 = 2.3 t
您实际上已经对所有可能的解决方案进行了参数化,并且因为您已经为它们分配了标识符 ln
,所以您不会 运行 陷入退化问题。即使任何 M[i][j]
是 t 的非线性函数,也没关系,因为系统本身是线性的,解完全由特征方程计算(在求解 lambda 时将保持 t 不变)。
我正在使用 GSL 库中的函数 gsl_eigen_nonsymm and/or gsl_eigen_symm 来查找 L x L 矩阵 M[i][j]
的特征值,它也是时间 t = 1,....,N
所以我有 M[i][j][t]
来获取每个 t 的特征值 我分配一个 L x L 矩阵 E[i][j] = M[i][j][t]
并为每个 t 对角化它。
问题是程序在一些迭代后以不同的顺序给出了特征值。例如 (L = 3) 如果在 t = 0
我在 t = 1
得到 eigen[t = 0] = {l1,l2,l3}(0)
我可能会得到 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");
}
}
我不认为你可以为所欲为。随着 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
中改变,那么你有一个解决方案,但我真的怀疑情况是否如此。
你的问题毫无意义。特征值对它们没有任何固有的顺序。在我看来,您想定义 M_t 的特征值类似于 L_1(M_t),..., L_n(M_t) 然后跟踪它们如何随时间变化。假设您的过程驱动 M_t 是连续的,那么您的特征值也是连续的。换句话说,当您对 M_t 进行小的更改时,它们不会发生显着变化。因此,如果您通过强制执行 L_1 < L_2... < L_n 来定义排序,那么此排序不会因 t 的微小变化而改变。当您有两个特征值交叉时,您需要决定如何分配更改。如果 "random fluctuations" 大于特征值之间的典型距离,那么这基本上是不可能的。
这是另一种跟踪特征向量的方法,可能会更好。为此,假设您的特征向量为 v_i,分量为 v_ij 。你首先要做的是 "normalize" 你的特征向量使得 v_i1 是非负的,即适当地翻转每个特征向量的符号。这将通过对 v_i1(每个特征向量的第一个分量)的排序来定义特征值的排序。这样您仍然可以跟踪相互交叉的特征值。但是,如果您的特征向量在第一个组件上交叉,您就有麻烦了。
根据文档,这些函数 return 无序:
https://www.gnu.org/software/gsl/manual/html_node/Real-Symmetric-Matrices.html
This function computes the eigenvalues of the real symmetric matrix A. Additional workspace of the appropriate size must be provided in w. The diagonal and lower triangular part of A are destroyed during the computation, but the strict upper triangular part is not referenced. The eigenvalues are stored in the vector eval and are unordered.
即使是 return 排序的函数,也可以通过简单的 ascending/descending 量级来实现:
https://www.gnu.org/software/gsl/manual/html_node/Sorting-Eigenvalues-and-Eigenvectors.html
This function simultaneously sorts the eigenvalues stored in the vector eval and the corresponding real eigenvectors stored in the columns of the matrix evec into ascending or descending order according to the value of the parameter sort_type as shown above.
如果您正在寻找特征值的时间演化,就像您一直在做的那样并求解时间相关的表示,例如:
lambda_1(t).....lambda_N(t)
对于简单的时间标量示例,
l1 = -1.3 t , l2 = -t , l3 = 2.3 t
您实际上已经对所有可能的解决方案进行了参数化,并且因为您已经为它们分配了标识符 ln
,所以您不会 运行 陷入退化问题。即使任何 M[i][j]
是 t 的非线性函数,也没关系,因为系统本身是线性的,解完全由特征方程计算(在求解 lambda 时将保持 t 不变)。