Tracking Eigenvalues Through a Crossing

Question

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?

Answer 1

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

Answer 2

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

Answer 3

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

Answer 4

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.

Here is your example:

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.

Answer 5

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

ListPlot[res]

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

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.

Answer 6

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

Answer 7

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.

Answer 8

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

---

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

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