# Tracking Eigenvalues Through a Crossing

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.

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 = {};
Do[
evalsN = Eigenvalues[H];
AppendTo[list, evalsN],
{B, 0, 0.1, 0.001}]


then I get a messed up plot like this:

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?

• 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"). – David G. Stork Feb 4 '18 at 21:00
• How large is large? – Michael E2 Feb 4 '18 at 21:03
• 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. – Jack S Feb 4 '18 at 21:16
• 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." – David G. Stork Feb 4 '18 at 22:03
• Do eigenvalues cross or do they bounce off in opposite directions? Eigenvalues don't have an identity to track. – alancalvitti Feb 6 '18 at 18:02

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

ClearAll[L];
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]]


• It's good to know there is some way to keep the integration from jumping paths, at least when the crossings are nontangential. – Daniel Lichtblau Feb 6 '18 at 19:23

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]


True

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


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

• Interesting approach, but this seems like it will break down for very large matrices. – Jack S Feb 6 '18 at 3:34
• @JackS, no idea how robust this approach is. – kglr Feb 6 '18 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,
lp,
lp[[o[[InversePermutation@n]]]]
]
]
];
sol = Quiet[
NDSolveValue[
{
\[FormalR]'[s] == lambdaprime[s, old[s]], old[s] == \[FormalR]'[s],
\[FormalR][xi] == lambda[xi], old[xi] == lambdap[xi]
},
\[FormalR],
{s, x0, x1}
],
NDSolveValue::pdord
]
]

EigenvaluePlot[matrix_, {x_, x0_, x1_}] := Module[{sol},
sol = ParametricEigenvalues[matrix, {x, x0, x1}];
Plot[
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.

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


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.

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},
x,
{B,-.1,.1},
Method->{"Projection","Invariants"->poly},
AccuracyGoal->10,
PrecisionGoal->10,
WorkingPrecision->20
];


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

Visualization:

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


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},
Do[
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 *)
Do[
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}]


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} *)


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 /@
Table[
Eigenvalues[H],
{B, 0, .1, .001}
];


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

projSortEVs[evs_] :=
Block[
{
new,
last = evs[[1]],
fd = evs[[2]] - evs[[1]]
},
Prepend[First@evs]@
Map[
Function[
new = Part[#, First /@ Nearest[(last + fd) -> "Index", #]];
fd = new - last;
last = new
],
Rest@evs
]
];
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 =
Table[
Eigenvalues[H],
{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.

• 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] – Daniel Lichtblau Mar 11 '18 at 14:56

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.

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 *)
Table[
λ = FindEigenvalue[H, 2 λ′ - λ′′]; (* linear extrapolate guess *)
λ′′ = λ′; λ′ = λ; (* update older and old results *)
{B, λ}
, {B, -0.1, 0.1, 0.001}]
, {λ0, λ0s}];

ListPlot[res]


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 SparseArrayStronglyConnectedComponents function as in this solution:

Hb[B_] = H[[#, #]] & /@ SparseArrayStronglyConnectedComponents[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}]]]


$$^*$$ 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 SparseArrayStronglyConnectedComponents, 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];


or

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 SparseArrayStronglyConnectedComponents block-diagonalization to H1`. A bit of experimentation may be required.