Suppose I have a matrix which depends on some parameter. I want to compute the eigenvalues as a function of this parameter, and then plot them. For example, I may have a matrix representing the Hamiltonian of a system in a magnetic field, and I want to plot the energies of the various states as a function of B-field. Sometimes the eigenvalues may cross each other, but I want to make sure the right eigenvalue stays associated with its own state. For example, consider the matrix

  H = {{1/2 + 21 B, 0, 0, 0, 0, 0}, {0, 1/2 + 7 B, 0, 0, 7 Sqrt[2] B, 
      0}, {0, 0, 1/2 - 7 B, 0, 0, 7 Sqrt[2] B}, {0, 0, 0, 1/2 - 21 B, 0, 
      0}, {0, 7 Sqrt[2] B, 0, 0, -1 + 14 B, 0}, {0, 0, 7 Sqrt[2] B, 0, 
      0, -1 - 14 B}}

If I evaluate evals = Eigenvalues[H] and then plot the result, I see a nice plot of the eigenvalues and the states follow properly through crossings, for example just below B = 0.05. Good crossing

As my matrix gets large, it becomes extremely slow to do the eigendecomposition analytically, so I'd rather do it numerically. However, in this case the eigenvalues don't properly track through crossings-- instead, they are sorted by value from largest to smallest in absolute value. For example, if I do the following

list = {};
 evalsN = Eigenvalues[H];
 AppendTo[list, evalsN],
 {B, 0, 0.1, 0.001}]

then I get a messed up plot like this:

Bad plot

This is for two reason: (1) because the ordering is all off, and (2) because the states don't track through crossing. I can fix problem (1) by ordering the eigenvalues by their magnitude, but that does NOT help them track through a crossing (e.g., see how the red and purple lines don't cross through one another around element 50). How do I do the second task?

  • $\begingroup$ Sort will work fine if the eigenvalues do not cross. You could also track the derivative of each eigenvalue with respect to your parameter and use that to associate eigenvalues at different values of the parameter (as this will remain nearly constant across a "crossing"). $\endgroup$ Commented Feb 4, 2018 at 21:00
  • 1
    $\begingroup$ The point of my question was that Sort will not work due to crossings. In reality, my matrix will be about 48x48 or greater in size, and have many (>10) crossings. $\endgroup$
    – Jack S
    Commented Feb 4, 2018 at 21:16
  • 2
    $\begingroup$ I would SortBy the product of the value of an eigenvalue and its derivative with respect to the parameter. That product will remain nearly constant across a "break." $\endgroup$ Commented Feb 4, 2018 at 22:03
  • 2
    $\begingroup$ Do eigenvalues cross or do they bounce off in opposite directions? Eigenvalues don't have an identity to track. $\endgroup$ Commented Feb 6, 2018 at 18:02
  • 1
    $\begingroup$ @alancalvitti I agree that there are caveats in what I wrote. However we're talking quantum physics here, diagonalization of Hamiltonian operators. There's rarely any singular or bouncy behavior there, and so maybe I took some heuristic shortcuts that aren't available in the mathematically general case. $\endgroup$
    – Roman
    Commented Jun 4, 2019 at 14:02

8 Answers 8


I'm not sure this will be robust enough for a 48 x 48 matrix.

The idea is common enough: Build an interpolation of a function by integrating its derivative with NDSolve, using a particular value of the function as an initial condition. If the derivative is computed as continuous, then the functions we get should track an eigenvalue trajectory through a crossing.

The eigenvalues satisfy known relationships with the coefficients of the characteristic polynomial of the matrix (the coefficients are elementary symmetric functions of the eigenvalues). We differentiate twice and project the solution onto the equations and their first derivatives. Projecting the first derivative of the solution onto the differentiated equations should make it unlikely that the solution will jump tracks, unless the tracks are tangent (which is an indeterminate case).

However, the OP's example, and perhaps all matrices, lead to an ODE that seems ill-conditioned near a point where an eigenvalue has a multiplicity greater than 1. Trying different methods and working precisions, sometimes the integration stopped at a crossing (too stiff) and at other times it jumped tracks at a crossing, even though that made the derivative discontinuous. Using the "FixedStep" method is the only way I could find to get the solution to stay on its tracks. It jumps across the crossing. The projection method keeps the solution accurate. I cannot extend the integration to B == 0, where the multiplicities of the eigenvalues are {4, 2}.

SetAttributes[L, NHoldAll];

(ode1 = D[
    eqn = SymmetricPolynomial[#, Array[L[#][B] &, Length@H]] & /@ Range@Length@H
        == Reverse@Most@CoefficientList[
          CharacteristicPolynomial[H, x], x] (-1)^Range@Length@H, {B, 1}];
 ode2 = D[ode1, B];
 ode = Solve[ode2, Array[L[#]''[B] &, Length@H]] /. Rule -> Equal;) // AbsoluteTiming

(*  {9.78918, Null}  *)

B0 = 1/10;
{sol} = NDSolve[{ode, 
     Array[L[#][B0] &, Length@H] == 
      Eigenvalues[H /. B -> B0], (ode1 /. B -> B0 /. 
       Thread[Array[L[#][B0] &, Length@H] -> Eigenvalues[H /. B -> B0]])}, 
    Array[L, Length@H], {B, 1/10000, 1/10}, 
    Method -> {"FixedStep", 
      Method -> {"Projection", "Invariants" -> {ode1, eqn}}}, 
    StartingStepSize -> 1/512, 
    WorkingPrecision -> 32]; // AbsoluteTiming

(*  {4.96584, Null}  *)

Plot[Evaluate[Array[L[#][B] &, Length@H] /. sol], 
 Evaluate@Flatten@{B, L[1]["Domain"] /. sol}, 
 PlotLegends -> Array[L, Length@H]]

Mathematica graphics

  • 7
    $\begingroup$ It's good to know there is some way to keep the integration from jumping paths, at least when the crossings are nontangential. $\endgroup$ Commented Feb 6, 2018 at 19:23

Here's a CharacteristicPolynomial-free approach based on the formula for eigenvalue sensitivity ${d\lambda \over dB}={\vec u {dH \over dB} \vec v \over <\vec u,\vec v>}$ where $\vec v$ is the corresponding right eigenvector and $\vec u$ is the corresponding left eigenvector (see e.g. eqn. (3.29) in Caswell 2019). We will start on an eigenvalue and track it using this formula in NDSolve. Because your matrix is symmetric, $\vec u$=$\vec v$.

First, a function to find the (right) eigenvector corresponding to a given eigenvalue, based on this answer by @mikado.

Eigenvector[mat_, λ_?NumericQ] := Module[{ev},
    ev = NullSpace[mat - λ IdentityMatrix[Length[mat]],
      Tolerance -> 10^tolp];
    If[Length[ev] == 1, Return[ev[[1]]]];
  , {tolp, -12, -2}]

Note that I gradually loosen the Tolerance until an answer is found, due to numerical errors in the eigenvalue. Also note that this gives a single eigenvector in the case of repeated eigenvalues. Surprisingly this doesn't cause any problems at the crossings.

Next, a function to compute the change in the eigenvalue with respect to the parameter (assuming symmetric matrix):

dλ[λ_?NumericQ, b_?NumericQ] := Module[
  {v = Eigenvector[H /. B -> b, λ]},
  v.(D[H, B] /. B -> b).v

No need to normalize thanks to NullSpace I think.

Finally, plug this function into NDSolve with initial values taken at a well-behaved parameter value.

B0 = 0.1; (* initial parameter value *)
λ0s = Eigenvalues[H /. B -> B0]; (* initial eigenvalues *)
  sol[i] = NDSolve[{λ'[B] == dλ[λ[B], B], λ[B0] == λ0s[[i]]},
    λ, {B, -0.1, 0.1}][[1]];
, {i, Length[λ0s]}];

Plot[Evaluate[λ[B] /. Table[sol[i], {i, Length[λ0s]}]], {B, -0.1, 0.1}]

Mathematica graphics

Frankly, I'm surprised this works so well. I thought the solution would easily get off track or break down at the intersections, but comparing values at the far side B==-0.1, looks OK.

N[Eigenvalues[H /. B -> -0.1]] // Sort
Table[λ[-0.1] /. sol[i], {i, 6}] // Sort
(* {-2.77986, -1.6, -0.267708, 0.179865, 1.86771, 2.6} *)
(* {-2.77986, -1.6, -0.267708, 0.179865, 1.86771, 2.6} *)
  • $\begingroup$ Please, clarify for me if this can be used with Eigensystem or numerically? I see you anticipate the result to be normalized without this, though...my system is composed of block matrices of the unperturbed system along the diagonals, with the perturbations elsewhere. The matrix is large (minimum of 20x20 of 2x2 blocks) & hermitian, so I obtain valuable complex & non-complex terms I need to keep. I do only plot Real parts, but I think this can work before I filter them? Looks like H would be the matrix with B dependence? Thanks in advance, I'm probably missing something super basic! $\endgroup$ Commented Aug 11, 2019 at 8:12
  • 1
    $\begingroup$ @CATrevillian Hard to say, why not open a new question with all the details? $\endgroup$
    – Chris K
    Commented Aug 11, 2019 at 8:23
  • $\begingroup$ can you have a new question be about a specific answer? That would change things, certainly, hmmm....I wouldn’t know how to generalize it, but maybe after the morning coffee brews. My problem is that I haven’t published yet, so I don’t want to give up the schrödinger’s cat, if you will...I’ll tag you if/when I do! I think this method will be very beneficial, though, thank you for this. Hopefully I can finagle it to work with non symmetric matrices!! $\endgroup$ Commented Aug 11, 2019 at 9:16
  • 1
    $\begingroup$ @CATrevillian For non-symmetric matrices, just calculate the left eigenvector $\vec u$ and put that in the formula for $d \lambda$. I'd just make the new question about your problem and see what people come up with. $\endgroup$
    – Chris K
    Commented Aug 11, 2019 at 9:29
  • $\begingroup$ Well, something truly strange is going on! I think it may work, but you don't define λ? Not sure if I lost something in translation, but for Hermitian matrices at the very least the left evec is the ConjugateTranspose of the right evec. I had to just take the ComplexConjugate to get past that error, but then it gives me Non-numeric errors for NDSolve. This presumably doesn't work with complex-valued problems? With time, I may figure it out...for now, just a headache has been produced! Many thanks, however, this should yield something good--eventually! $\endgroup$ Commented Aug 11, 2019 at 17:53

We can avoid the issue working with the roots of CharacteristicPolynomial of H:

cproots = x /. Solve[CharacteristicPolynomial[H, x] == 0, x];
Reduce[cproots == Eigenvalues[H], Reals]


list1 = Transpose @ Table[cproots, {B, 0, .1, .001}];
ListLinePlot[list1, BaseStyle -> Thick, DataRange -> {0, .1}]

enter image description here

Using cproots in place of Eigenvalues[H] in OP's code gives the same result:

list = {};
Do[evalsN = cproots; AppendTo[list, evalsN], {B, 0, 0.1, 0.001}]
ListLinePlot[Transpose@list, BaseStyle -> Thick, DataRange -> {0, .1}]

same picture

  • $\begingroup$ Interesting approach, but this seems like it will break down for very large matrices. $\endgroup$
    – Jack S
    Commented Feb 6, 2018 at 3:34
  • $\begingroup$ @JackS, no idea how robust this approach is. $\endgroup$
    – kglr
    Commented Feb 6, 2018 at 19:51

Here's an alternate approach also using NDSolve. The basic idea is that at a crossing, the order of the eigenvalues may change. This means that the order of the eigenvalue derivatives will also change. If we assume that at a crossing the eigenvalue derivatives are different, then we can try to correct the order by looking at the Ordering of the derivatives. Here is a function that does this:

ParametricEigenvalues[matrix_, {x:Except[_List], x0_, x1_}] := 
ParametricEigenvalues[matrix, {{x, x0}, x0, x1}]

ParametricEigenvalues[matrix_, {{x_, xi_}, x0_, x1_}] := Module[
    {lambda, lambdap, lambdaprime, sol, old},
    lambda = Function[x, Evaluate[Eigenvalues[matrix]]];
    lambdap = lambda';
    lambdaprime[s_, old_List] := With[{lp = lambdap[s]},
        With[{o = Ordering[old], n = Ordering[lp]},
            If[n === o,
    sol = Quiet[
                \[FormalR]'[s] == lambdaprime[s, old[s]], old[s] == \[FormalR]'[s],
                \[FormalR][xi] == lambda[xi], old[xi] == lambdap[xi]
            {s, x0, x1}

EigenvaluePlot[matrix_, {x_, x0_, x1_}] := Module[{sol},
    sol = ParametricEigenvalues[matrix, {x, x0, x1}];
        Evaluate @ Table[
            Indexed[sol[t], i],
            {i, Last @ Dimensions @ sol["ValuesOnGrid"]}
        {t, x0, x1}

Note that I include an auxiliary non-differential equation old[s] == \[FormalR]'[s] so that I can include the old derivative in my lambdaprime function.

Here is your example:

EigenvaluePlot[H, {B, -.1, .1}]

enter image description here

Note how the eigenvalues are correctly tracked at both $x=0$ and $x=0.047`$.

However, I think my idea of using derivative Ordering is flawed, so I would welcome hearing about examples where it doesn't work. Probably I should also be looking at eigenvalue ordering as well.


Here's another CharacteristicPolynomial-free approach, this one using Eigenvalues[Method -> "Arnoldi"] to find a single eigenvalue, then tracking it over parameter values using linear extrapolation (or NDSolve, see addendum).

First, wrap the Arnoldi eigenvalue solver in a function to find one eigenvalue given an initial guess.

FindEigenvalue[mat_, ic_?NumericQ] := Check[
  Eigenvalues[mat, 1, Method -> {"Arnoldi", Shift -> ic}][[1]], 
  ic, Eigenvalues::ssing];

Next, generate a list of initial solutions at a nice parameter value.

λ0s = Eigenvalues[H /. B -> -0.1]

Finally, do an outer Table over those initial solutions and an inner Table over the range of parameter values.

res = Table[
  λ′′ = λ′ = λ0; (* initialize older and old results *)
    λ = FindEigenvalue[H, 2 λ′ - λ′′]; (* linear extrapolate guess *)
    λ′′ = λ′; λ′ = λ; (* update older and old results *)
    {B, λ}
  , {B, -0.1, 0.1, 0.001}]
, {λ0, λ0s}];


Mathematica graphics

Addendum June 2023:

Instead of the final Table, another way to vary the parameter is with NDSolve (DAE mode):

sol = NDSolve[Flatten@Join[{x'[B] == 0},
  Table[{λ[n][B] == FindEigenvalue[H, λ[n][B]], λ[n][-0.1] == λ0s[[n]]}, {n, 6}]], 
  Table[λ[n], {n, 6}], {B, -0.1, 0.1}][[1]];

Plot[Evaluate[Table[λ[n][B] /. sol, {n, 6}]], {B, -0.1, 0.1}]

enter image description here

where x is a dummy variable to convince NDSolve to help us.

One note (not relevant to this problem): this fails if the eigenvalues are complex. You can work around that with ReIm, but run into another problem of stiffness if the eigenvalues change from real to complex conjugates.


If it is feasible to construct the characteristic polynomial of the matrix, then here is another NDSolve approach that works, although one might need to adjust some of the NDSolve options.

poly = CharacteristicPolynomial[H, x[B]]
init = Flatten @ Values @ Solve[poly == 0 /. B->1/10, x[1/10]]
sol = NDSolveValue[
    {D[poly, B] == 0, x[1/10] == init},

1/4 (1/2 - 21 B - x[B]) (1/2 + 21 B - x[B]) (1 - 2 x[B] - 3 x[B]^2 - 1764 B^2 x[B]^2 + 4 x[B]^3 + 4 x[B]^4)

{-(8/5), 13/5, 1/10 (8 - Sqrt[114]), 1/10 (8 + Sqrt[114]), 1/10 (-13 - Sqrt[219]), 1/10 (-13 + Sqrt[219])}


Plot[Evaluate @ Table[Indexed[sol[x], i], {i, 6}], {x, -.1, .1}]

enter image description here


If we block-diagonalize the matrix H, we can calculate the eigenvalues of the blocks separately. These eigenvalues often don't have crossings.

In the absence of additional info, we can do a simple$^*$ automated block-diagonalization using the SparseArray`StronglyConnectedComponents function as in this solution:

Hb[B_] = H[[#, #]] & /@ SparseArray`StronglyConnectedComponents[H]

(*  { {{1/2 + 21 B}},
      {{1/2 + 7 B, 7 Sqrt[2] B}, {7 Sqrt[2] B, -1 + 14 B}},
      {{1/2 - 7 B, 7 Sqrt[2] B}, {7 Sqrt[2] B, -1 - 14 B}},
      {{1/2 - 21 B}} }  *)

Calculate the eigenvalues of each block separately, and join:

ev[B_?NumericQ] := Join @@ Sort@*Eigenvalues /@ Hb[N[B]]

There are no crossings now:

ListLinePlot[Transpose[Table[Thread[{B, ev[B]}], {B, 0, 0.1, 0.001}]]]

enter image description here

$^*$ This is not a proper block-diagonalization, as some of the resulting blocks can potentially be further block-subdivided by a judicious unitary transformation. Inspiration for such unitary transformations can come from many sources, for example physical insight on conserved quantities (total angular momentum, etc.).

If your matrix has no discernible block structure for SparseArray`StronglyConnectedComponents, you may still be able to pre-diagonalize it by using a unitary transformation related to its eigenvectors, for example

U = Normalize /@ Eigenvectors[H /. B -> 0];


U = Normalize /@ Eigenvectors[D[H, B] /. B -> 0];

or similiar. Then, you can transform the matrix H to

H1 = Chop[Expand[U.H.ConjugateTranspose[U]]];

and apply the SparseArray`StronglyConnectedComponents block-diagonalization to H1. A bit of experimentation may be required.

  • $\begingroup$ When you say this is not proper block-diagonalization, and give the requisite advisories, do you refer to the possibility that it would not work at all? Or just not complete enough to be fully block-diagonalized? $\endgroup$ Commented Aug 18, 2019 at 22:15
  • 1
    $\begingroup$ @CATrevillian sometimes it won't work at all. The strongly connected components refer to structural blocks in the matrix, which can be identified only if the matrix can be block-diagonalized simply by reordering the rows/columns. This is very different from a mathematical block-diagonalization, which uses a much larger group of transformations including orthogonal and unitary transformations that disentangle structurally connected blocks. $\endgroup$
    – Roman
    Commented Aug 19, 2019 at 9:16
  • $\begingroup$ This makes sense, so presuming one has assembled the matrix using iterations of two sets of blocks, one for the diagonal components, and one for the not-diagonal components, should one expect to see any beneficial reorganization occur? Interaction between said blocks then will not occur, only pure structural reorganization? OP’s example used in this answer seems to suggest elimination of some rows/columns occur, but this is just due to their being (Real) zeroes in the structure? Thank you for your experience and knowledge with this stuff, Roman! $\endgroup$ Commented Aug 19, 2019 at 10:21
  • 1
    $\begingroup$ Yes it's the zeros in the structure that makes the structural component analysis possible. If your matrix is dense, then any structure is not visible to this analysis, and you need to resort to more mathematical tools. $\endgroup$
    – Roman
    Commented Aug 19, 2019 at 12:53
  • $\begingroup$ +1 for you, -1 for me haha—I’ve got a Dense Hermitian Matrix D; $\endgroup$ Commented Aug 19, 2019 at 13:20

Here's a way that sort-of works, based on the projecting the new value from the first-difference. It requires that we first sort our original list of EVs:

evs =
  Sort /@
    {B, 0, .1, .001}

Then we do a moving ordering based on the projected value from the previous difference:

projSortEVs[evs_] :=
    last = evs[[1]],
    fd = evs[[2]] - evs[[1]]
      new = Part[#, First /@ Nearest[(last + fd) -> "Index", #]];
      fd = new - last;
      last = new
projSortEVs[evs] // Transpose // ListLinePlot[#, DataRange -> {0, .1}] &


Note that this will still be slow since we build a NearestFunction at each step. Unfortunately I don't have a better way to do that.

It also fails miserably if we start from unsorted data:

evsUS =
   {B, 0, .1, .001}
projSortEVs[evsUS] // Transpose // 
 ListLinePlot[#, DataRange -> {0, .1}] &


Why this first starts to go bad for such linear-looking data I'm not sure (also note that we don't lose the ordering until like the 16th point). Once it goes bad obviously it can't be fixed since the next step depends on the previous.

In general, it might only be able to correct pretty obvious errors.

  • 2
    $\begingroup$ Might be able to speed it by only getting careful near crossing points. If cp is the characteristic polynomial in the variable t, these can be found as NSolve[Discriminant[cp, t] == 0, B] $\endgroup$ Commented Mar 11, 2018 at 14:56

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.