Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I've been working on this research problem for months with some success by no global/proper solution to the problem.

So I have square, Hermitian matrices of various sizes ($8,20,38$ dimensions) with functional dependence on two variables, call them x and y. I am interested in plotting all the eigenvalues as functions of the $2$ variables x and y.

Let's call M the $8 \times 8$ hermitian matrix.

Plot3D[ Eigenvalues[M], {x,-1,1}, {y,-1,1}]

gives $8$ sheets corresponding to the eigenvalues of M as functions of x and y.

However, there will be many "jumps" between the sheets which totally obscure the plots. Basically, I think that during the plotting of the eigenvalues, for certain regions of the parameter space, Mathematica organizes the $8$ values differently, so the different sheets swap places, so we get tons of random garbage vertical bits between sheets.

I have spent months trying to resolve this in different ways, mostly involving methods to force Mathematica into a consistent sorting method. Here are some examples of failed methods:

Plot3D[ Sort[Eigenvalues[M]], {x,-1,1}, {y,-1,1}]

Plot3D[ Sort[Eigenvalues[M], Abs[#1] <= Abs[#2]&], {x,-1,1}, {y,-1,1}]

Block[{set=Ordering[ Eigenvalues[ M/.x->0 /. y->0], All, Abs[#1] <= Abs[#2]&]},
  Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]

Block[{set=Ordering[ Ordering[Eigenvalues[ M/.x->0 /.y->0], All, 
        Abs[#1] <= Abs[#2]&]]},
  Plot3D[ Eigenvalues[M][[set]], {x,-1,1}, {y,-1,1}]]

Block[{set = Ordering[Eigenvalues[ M /.x->0 /.y->0], All, Abs[#1] <= Abs[#2]&]},
  Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]

Block[{set = Ordering[ Ordering[ Eigenvalues[ M /.x->0 /.y->0], All,
       Abs[#1] <= Abs[#2]&]]},
  Plot3D[ Sort[ Eigenvalues[M][[set]]], {x,-1,1}, {y,-1,1}]]

as well as other permutations of similar ideas and more. Some of these methods have worked somewhat, reducing the jumping or working in certain regimes of the parameter domain. But nothing works consistently.

Please help me understand the proper solution to this problem!

share|improve this question
do you have an actual example? –  Nasser Jan 3 at 17:14
Not good: i.imgur.com/aSjHES5.png?1 –  Steve Jan 3 at 17:21
good: i.imgur.com/vWeSYGc.png?1 –  Steve Jan 3 at 17:22
Example: m = SparseArray[{i_, j_} -> Sin[i j x y], {3, 3}]; Plot3D[ Eigenvalues[m][[1]], {x, 0, 1}, {y, 0, 1}, PlotRange -> All] –  belisarius Jan 3 at 18:13
Try generating a Table, then working on an algorithm to sort the table, finally using ListPlot3D.. –  george2079 Jan 3 at 18:20

2 Answers 2

note see bottom of answer for a required fix for version 10+

A 2D example of what i suggested in my comment:

Based on Belisarius example..we get three lines obviously not connected properly:

 m = SparseArray[{i_, j_} -> Sin[i j 9 /10 y], {3, 3}];
 alle = Table[Eigenvalues[m], {y, 0, 1, .1}];
 original = Show[
 MapIndexed[ListPlot[Flatten[Take[alle, All, {#}] , 2],
      Joined -> True, PlotRange -> All, 
      PlotStyle -> Hue[First@#2/3]] &, Range[3]]]

enter image description here

Now construct a sort of surface curvature penalty function and minimize by selecting from each of the possibilities..

v = Table[V[i], {i, Length[alle]}];
cons  = Table[1 <= V[i] <= 3, {i, Length[alle]}];
SetOptions[NMinimize, MaxIterations -> 500]                    
v1 = v /. Last@
     Minimize[ {Total[((#[[1]] + #[[3]] - 2 #[[2]])^2 & /@ 
        Partition[MapIndexed[#[[v[[First@#2]]]] &, alle] , 3, 
           1, {1, 3}, {}])], cons}, v, Integers]

Nicely finds a smooth line..

enter image description here

   ListPlot[MapIndexed[#[[v1[[First@#2]]]] &, alle], Joined -> True, 
        PlotStyle -> {Thick, Black, Dashed}]]

Repeat the process adding to the constraints like this:

cons  = Join[cons, Table[V[i] != v1[[i]], {i, Length[alle]}] ]

enter image description here

Works well even if the input ordering is a random mess:

enter image description here

Quite likely horribly slow in 3d...though you might tackle it as a series of 2d sorts and join them together.

See here for that SetOptions[NMinimise] trick..


Edit - a more compact version, reordering the points instead of adding more constraints at each step:

 v = Table[V[i], {i, Length[alle]}];
 SetOptions[NMinimize, MaxIterations -> 500];
      alle = MapThread[ Prepend[  Drop[ #2 , {#1}]  , #2[[#1]] ]  &,
           { (v /. Last@Minimize[{Total[((#[[1]] + #[[3]] - 2 #[[2]])^2 & /@
        Partition[MapIndexed[#[[v[[First@#2]]]] &, alle], 3, 1, {1, 3}, {}])],
              Table[k <= V[i] <= 3, {i, Length[alle]}]}, v, Integers]) , 
        alle}] , {k, 2}];
    ListPlot[Flatten@Take[alle, All, {#}], Joined -> True, 
    PlotStyle -> {Thick, Black, Dashed}] & /@ Range[3]]

v10+ fix

for mathematica versions with the Indexed function, replace :

 MapIndexed[#[[v[[First@#2]]]] &, alle]


 Indexed @@@ Transpose[{alle, v}]

(and setting Off[Part::pspec]; is no longer needed )

share|improve this answer
Haven't see it before. Very nice –  belisarius Mar 3 at 20:38
Sounds great. However, with mathematica v10, I get the following error: Part::pkspec1: The expression V[1] cannot be used as a part specification. >> Is there something wrong? after the instuction v1 = v /. Last@ Minimize[ ..... –  altroware Aug 13 at 13:52
i guess the warning i surpressed got elevated to an error. i dont have 10 to work on it though –  george2079 Aug 13 at 14:19
@altroware fixed that issue i think - verified only with 9 which it turns out has the Indexed function (undocumented). –  george2079 Aug 19 at 17:57
up vote 2 down vote accepted

Thanks to everyone for the help. Here's my working solution:

First, the spirit of the underlying problem and the nature of the ultimate solution. Simply sorting by the function values (executed with the default one-argument Sort command) is enough to separate the functions and prevent the jumping. However, it does not give the user control over the functions and it is difficult to do anything further (e.g. remove several of the functions from the plot, apply colors in a predictable way, etc.). If you sort instead by the two argument Sort with a pure absolute value function as


there is a pairwise degeneracy problem when two functions cross planes equidistant above and below the zero plane, and you get tons of pairwise swapping at those places and a garbage plot. The spirit of this attempt, however, would be to gain control over the smallest absolute value function bands, for example.

The solution that follows uses a different technique. We instead use a reference point (the origin) to compute a single list of eigenvalues, perform simple operations there to compute which bands in the final plot will correspond to meeting a desired criteria (such as the smallest n bands in terms of absolute value) and then plot the functions sorted by a simple Sort[] with these indices

"lownumber" is the desired number of plotted eigenvalues, the smallest lownumber bands in terms of absolute value. lowb and lowt are the lower and upper bound indices for the plotted functions

     Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All, 
     Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];  

     Sort[Eigenvalues[M/. x -> 0 /. y -> 0]], All, 
     Abs[#1] <= Abs[#2] &], -(Dimensions[M][[1]] - lownumber)]]];

After this, there are two possible solutions. I have tested them both to be of the same speed to within about 3%, for matrices up to size 38. For much larger matrices, I am not sure, one might be preferable. Here they are:

 Table[Sort[Eigenvalues[hamiltonian]][[s]], {s, lowb, lowt}] 


 Table[RankedMin[Eigenvalues[hamiltonian],s], {s, lowb, lowt}] 
share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.