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

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

• 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! – CA Trevillian Aug 11 at 8:12
• @CATrevillian Hard to say, why not open a new question with all the details? – Chris K Aug 11 at 8:23
• 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!! – CA Trevillian Aug 11 at 9:16
• @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. – Chris K Aug 11 at 9:29
• 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! – CA Trevillian Aug 11 at 17:53

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


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

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.

• 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? – CA Trevillian Aug 18 at 22:15
• @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. – Roman Aug 19 at 9:16
• 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! – CA Trevillian Aug 19 at 10:21
• 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. – Roman Aug 19 at 12:53
• +1 for you, -1 for me haha—I’ve got a Dense Hermitian Matrix D; – CA Trevillian Aug 19 at 13:20