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!

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

3 Answers 3


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

see here for an alternate approach

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 )

  • $\begingroup$ Haven't see it before. Very nice $\endgroup$ Commented Mar 3, 2014 at 20:38
  • $\begingroup$ 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[ ..... $\endgroup$
    – altroware
    Commented Aug 13, 2014 at 13:52
  • $\begingroup$ i guess the warning i surpressed got elevated to an error. i dont have 10 to work on it though $\endgroup$
    – george2079
    Commented Aug 13, 2014 at 14:19
  • $\begingroup$ @altroware fixed that issue i think - verified only with 9 which it turns out has the Indexed function (undocumented). $\endgroup$
    – george2079
    Commented Aug 19, 2014 at 17:57
  • $\begingroup$ @george2079 would you care to help me out with my question where I try to use your above answer (the questions are strongly related): mathematica.stackexchange.com/questions/83906/… $\endgroup$
    – NOhs
    Commented May 27, 2015 at 8:01

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

I know the question was answered long ago but I ran into the same problem and was going to implement the solution above when I may have stumbled on a much simpler solution. I would run the following code for a 16x16 matrix (which is block diagonal with the largest blocks being 4x4)

Plot[Eigenvalues[FullHamiltonian], {x, -1, 1}]

and would obtain the horrible graph

A lot of jumps between branches.

I then realised that when I tried to plot a single eigenvalue I could use the horrible Root["polynomial"&,1] expression for a given eigenvalue to obtain just one of these lines in a smooth way. I therefore assumed it might be possible to solve it using

Plot[EV, {x,-1,1}]

And in fact this gives me the following plot

The jumps are all gone.

I repeated the experiment using two variables (on a 4x4 matrix) and it still seems to work. Jumps in 3D.Smooth in 3D.

My understanding is that the Root function has a well defined ordering whereas the Eigensystem function doesn't if all coefficients are numerical. When we do the diagonalisation analytically first, we are guaranteed to be using Root with the same argument every time. This should probably be slower than doing the diagonalisation of the fully numerical matrix but has the advantage of making the code much more legible.

As I have only tested it on an effectively 4x4 matrix it would be interesting to see what happens for larger matrices but I cannot see any obvious reason why it would fail.

EDIT: I tested on a 5x5 (first one without analytic solutions) and it still seems to work, however computing the eigenvectors analytically is very slow. There might be a way to combine the analytic eigenvalues, which are fast, with a fully numeric finding of the eigenvectors to combine speed and code sobriety.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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