# Sorting eigenvectors according to its projection

The problem

I'm trying to calculate the eigenvalues and eigenvectors of a matrix that depends on a parameter x. As x changes, I get a series of eigenvalues and eigenvectors corresponding to different x. Then I need to sort the eigenvalues and eigenvectors so that I get a continuous change. I can make it work, but the code runs too slow. Here are the details:

This is the function that generates the matrix that depends on x

H[x_,nBlock_:21] := Module[{tridiag, diagL, diag},
tridiag = {{0, x/2}, {x/2, 0}};
diag = {{-1/2, 0}, {0, 1/2}};
diagL = IdentityMatrix[{2, 2}];
SparseArray[{Band[{1, 1}] ->
Table[diag + diagL*(i - nBlock/2), {i, 1/2, nBlock - 1/2}],
Band[{3, 1}] -> Table[tridiag, {nBlock - 1}],
Band[{1, 3}] -> Table[tridiag, {nBlock - 1}]}]
]


For a specific x, we calculate the eigenvalues and eigenvectors, and sort them according to the eigenvalue. Then we do this for a series of x values.

{valls, vecls} =
Transpose@Table[
Transpose@SortBy[Transpose[Eigensystem[H[x]]], First], {x, 0.1, 5, 0.1}];


Here is what the 19th to 24th eigenvalues look like, as a function of x. We can see that there are some "avoided" crossings in the eigenvalues (notice how the colors change near those crossings).

ListPlot[Transpose[valls][[19 ;; 24]], Joined -> True]


And the goal is to fix that by reordering the eigenvalues and eigenvectors in each x so that the eigenvalues look like this (notice the colors of the curve)

ListPlot[Transpose[valSortls][[19 ;; 24]], Joined -> True]


The algorithm for reorderring is like this: Say we have eigenvectors {v1,v2,v3,...} corresponding to the matrix H[x], and {u1,u2,u3...} corresponding to H[x+dx], and we want to reorderring {u1,u2,u3...} based on {v1,v2,v3,...}. Then the correct vector at the first place should be the one in {u1,u2,u3...} that has the largest projection on v1. For example, if u1 is the correct vector, then it should satisfy Abs[Conjugate[v1].u1]>Abs[Conjugate[v1].u2] and Abs[Conjugate[v1].u1]>Abs[Conjugate[v1].u3]. The same is true for the other vectors in {u1,u2,u3...}. Moreover, I know that the crossing only happens between the neighbors.

My implementation

I tried a rudimentary implementation:

a function sort two vector lists, and return the ordering index:

sortVecs[vecls1_, vecls2_] := Module[{vecSort, lth, veclth, neib},
lth = Length[vecls1];
Table[
neib = Select[{n - 1, n, n + 1}, lth >= # >= 1 &];
First@Last[
SortBy[Transpose[{neib, vecls2[[neib]]}],
Abs[Conjugate[vecls1[[n]]].#[[2]]] &]]
, {n, 1, lth}]
]


and the function that sorts both the eigenvectors and eigenvalues:

sortEigensystem[{valls_, vecls_}] :=
Module[{lth = Length[valls], vallsSorted, veclsSorted, vec, val, ordering},
vallsSorted = {};
veclsSorted = {};
vec = vecls[[1]];(*list of vectors being sorted*)
val = valls[[1]];(*the corresponding values being sorted*)
veclsSorted = Append[veclsSorted, vec];
vallsSorted = Append[vallsSorted, val];

Table[

ordering = sortVecs[vec, vecls[[n + 1]]];
vec = vecls[[n + 1]][[ordering]];
val = valls[[n + 1]][[ordering]];

veclsSorted = Append[veclsSorted, vec];(*append the sorted vectors*)
vallsSorted = Append[vallsSorted, val];(*append the sorted values*)

, {n, 1, lth - 1}];

{vallsSorted, veclsSorted}

]


Here is a test

{valSortls, vecSortls} = sortEigensystem[{valls, vecls}]; // AbsoluteTiming
(* {0.030064, Null} *)


Question

How can I make sortEigensystem run faster? For my real situation, I need to sort about 1000 400X400 matrixes.

{valls, vecls} =
Transpose@
Table[Transpose@SortBy[Transpose[Eigensystem[H[x,201]]], First], {x,
0.1, 10, 0.01}];
{valSortls, vecSortls} =
sortEigensystem[{valls, vecls}]; // AbsoluteTiming
(* {27.7415, Null} *)


I would like to reduce the time to less than 10 seconds if that is possible.

• You may try this solution mathematica.stackexchange.com/a/39754/193 .Not sure if the question is an exact dup – Dr. belisarius Oct 1 '15 at 6:59
• @xslittlegrass - in your example, it takes only 30 milliseconds, which is not long. I assume you need to do this on either much matrices (larger than 42 by 42) or for many more x-values than 50. What are you looking to do with this in the end, that this is the choke-point timewise? – Jason B. Oct 1 '15 at 9:09
• Possible duplicate of How get eigenvectors without phase jump? – Daniel Lichtblau Oct 1 '15 at 14:47
• Appears to be a duplicate of 1, 2, and 3 (and maybe others). – Daniel Lichtblau Oct 1 '15 at 14:48
• @DanielLichtblau I think this is different than How get eigenvectors without phase jump since there we are dealing with the phase of the eigenvectors, whereas we are dealing with the orders of the eigenvectors here. But I agree the question is the same as in 3. – xslittlegrass Oct 1 '15 at 14:54

I just tried to clean up the code a bit. With respect to sortVecs, I used Ordering which is really what you were going for. You were wasting a little bit of time by taking the conjugate of one of the vectors, when the eigenvectors in vecls are all real-valued. Should you move to a different form for H that gives complex eigenvectors, just uncomment the relevant part below.

In sortEigensystem you were using x=Append[x,newval], which is equivalent to AppendTo[x,newval] which should be avoided as it is a time waste - (I think because at each step it needs to find the end of the list and then tack on a new bit). Use Reap and Sow instead (or even Table if possible, but in this case Table is slower than Reap and Sow).

Those changes made it about 2 to 2.5 times faster, depending on the size of the system, but I was able to make it much faster by compiling the sortVecs function.

Here are the new versions of sortVecs and sortEigensystem:

sortVecsComp = Compile[{{vecls1, _Real, 2}, {vecls2, _Real, 2}},
Module[{vecSort, lth, neib}, lth = Length[vecls1];
Table[neib = Select[{n - 1, n, n + 1}, lth >= # >= 1 &];
neib[[Last@Ordering[Abs[((*Conjugate@*)vecls1[[n]]).#] & /@
vecls2[[neib]]]]], {n, 1, lth}]]
];

sortEigensystem2[{valls_, vecls_}] :=
Module[{lth = Length[valls], vallsSorted, veclsSorted, vec, val,
ordering,va,ve},
vec = vecls[[1]];(*list of vectors being sorted*)

val = valls[[1]];(*the corresponding values being sorted*)
Reap[
Do[
ordering = sortVecsComp[vec, vecls[[n]]];
vec = vecls[[n]][[ordering]];
val = valls[[n]][[ordering]];
Sow[val, va];
Sow[vec, ve];
, {n, 1, lth}]][[2]]
];


Here is an example on $400\times 400$ matrices,

(* I have to use Quiet here because it shouts errors about the fact that it's
easier to run Eigensystem on normal arrays than on SparseArrays of this size *)
Quiet[{valls, vecls} =
Transpose@Table[Transpose@SortBy[Transpose[Eigensystem[H[x, 201]]], First], {x, 0.1, 10, 0.01}];]

{valSortls, vecSortls} =
sortEigensystem[{valls, vecls}]; // AbsoluteTiming
{valSortls2, vecSortls2} =
sortEigensystem2[{valls, vecls}]; // AbsoluteTiming
{valSortls2, vecSortls2} == {valSortls, vecSortls}

(* {25.1059, Null} *)

(* {1.83797, Null}  *)

(* True   *)


Hope that helps.

• Thanks for the answer. I didn't get such a large speedup as you did. I got 2.5X speedup, which is already very good. – xslittlegrass Oct 5 '15 at 16:36
• Is that before or after I edited and made changes to 'sortVecs'? – Jason B. Oct 5 '15 at 17:20
• After your edit on Oct 2. – xslittlegrass Oct 5 '15 at 18:03
• Just tried it again, and the dramatic results I had before weren't reproducible. Now I would only get a speedup of 2 to 3 times faster. But I couldn't let that be, so I tried to compile the sortVecs function, and it's loads faster now. – Jason B. Oct 6 '15 at 7:26
• OK, that's impressive. I was too focused onto the form of the programing and forgot the Compile trick. Thanks! – xslittlegrass Oct 6 '15 at 15:27