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.