GSL eigenvalues order

Question

I am using the functions gsl_eigen_nonsymm and/or gsl_eigen_symm from the GSL library to find the eigenvalues of an L x L matrix M[i][j] which is also a function of time t = 1,....,N so i have M[i][j][t] to get the eigenvalues for every t i allocate an L x L matrix E[i][j] = M[i][j][t] and diagonalize it for every t.

The problem is that the program gives the eigenvalues at different order after some iteration. For example (L = 3) if at t = 0 i get eigen[t = 0] = {l1,l2,l3}(0) at t = 1 i may get eigen[t = 1] = {l3,l2,l1}(1) while i need to always have {l1,l2,l3}(t) To be more concrete: consider the matrix M (t) ) = {{0,t,t},{t,0,2t},{t,t,0}} the eigenvalues will always be (approximatevly) l1 = -1.3 t , l2 = -t , l3 = 2.3 t When i tried to diagonalize it (with the code below) i got several times a swap in the result of the eigenvalues. Is there a way to prevent it? I can't just sort them by magnitude i need them to be always in the same order (whatever it is) a priori. (the code below is just an example to elucidate my problem)

EDIT: I can't just sort them because a priori I don't know their value nor if they reliably have a structure like l1<l2<l3 at every time due to statistical fluctuations, that's why i wanted to know if there is a way to make the algorithm behave always in the same way so that the order of the eigenvalues is always the same or if there is some trick to make it happen.

Just to be clearer I'll try to re-describe the toy problem I presented here. We have a matrix that depends on time, I, maybe naively, expected to just get lambda_1(t).....lambda_N(t), instead what I see is that the algorithm often swaps the eigenvalues at different times, so if at t = 1 I've got ( lambda_1,lambda_2,lambda_3 )(1) at time t = 2 (lambda_2,lambda_1,lambda_3)(2) so if for instance I wanted to see how lambda_1 evolves in time I can't because the algorithm mixes the eigenvalues at different times. The program below is just an analytical-toy example of my problem: The eigenvalues of the matrix below are l1 = -1.3 t , l2 = -t , l3 = 2.3 t but the program may give me as an output(-1.3,-1,2.3)(1), (-2,-2.6,4.6)(2), etc As previously stated, I was wondering then if there is a way to make the program order the eigenvalues always in the same way despite of their actual numerical value, so that i always get the (l1,l2,l3) combination. I hope it is more clear now, please tell me if it is not.

#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

Your question makes zero sense. Eigenvalues do not have any inherent order to them. It sounds to me like you want to define eigenvalues of M_t something akin to L_1(M_t),..., L_n(M_t) and then track how they change in time. Assuming your process driving M_t is continuous, then so will your eigenvalues be. In other words they will not significantly change when you make small changes to M_t. So if you define an ordering by enforcing L_1 < L_2... < L_n, then this ordering will not change for small changes in t. When you have two eigenvalues cross, you'll need to make a decision about how to assign changes. If you have "random fluctuations" which are larger than the typical distance between your eigenvalues, then this becomes essentially impossible.

Here's another way of tracking eigenvectors, which might prove better. To do this, suppose that your eigenvectors are v_i, with components v_ij . What you do is first "normalize" your eigenvectors such that v_i1 is nonnegative, i.e. just flip the sign of each eigenvector appropriately. This will define an ordering on your eigenvalues through an ordering on v_i1, the first component of each eigenvector. This way you can still keep track of eigenvalues that cross over each other. However if your eigenvectors cross on the first component, you're in trouble.

Answer 2

I don't think you can do what you want. As t changes the output changes.

My original answer mentioned ordering on the pointers but looking at the data structure it won't help. When the eigenvalues have been computed the values are stored in E. You can see them as follows.

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]);

The following code is how the the data in the eigenvector is layed out...

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

If, and you can check this when you see the order change, the order of the data in mdata changes AND the order in data changes then the algorithm does not have a fixed order ie you cannot do what you're asking it to do. If the order does not change in mdata and it changes in data then you have a solution but I really doubt that will be the case.

Answer 3

According to the docs, those functions return unordered:

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.

Even the functions that return ordered results, do so by simple ascending/descending magnitude:

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.

If you're looking for the time evolution of the eigenvalues, just do like you have been doing and solve for the time-dependent representations, e.g.:

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

For your simple time-as-scalar example,

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

You literally have a parameterization of all possible solutions and because you've assigned them identifiers ln you don't run into the issue of degeneracy. Even if any M[i][j] are nonlinear functions of t, it shouldn't matter because the system itself is linear and solutions are computed purely by the characteristic equation (which will hold t constant while solving for lambda).