Sorting eigenvectors according to its projection

Question

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.

Answer 1

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.