# Different result from utilizing of eigenvalues and eigenvectors commands

I have a matrix as a function of a parameter like

d[\[Delta]_] := {{(
6.20074189109658 -
0.5773502691896257 \[Delta])/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6),
0., -(0.19595917942265423/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6)),
0.27712812921102037/(-(3 + 2 \[Delta])^2 (-3 + 4 \[Delta]))^(
1/6), 0., 0.}, {0., (
6.339305955702091 +
1.1547005383792515 \[Delta])/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6), 0., 0.,
0.27712812921102037/(-(3 + 2 \[Delta])^2 (-3 + 4 \[Delta]))^(
1/6), 0.}, {-(
0.19595917942265423/(-(3 + 2 \[Delta])^2 (-3 + 4 \[Delta]))^(
1/6)), 0., (
1.1547005383792515 (5.25 + \[Delta]))/(-(3 +
2 \[Delta])^2 (-3 + 4 \[Delta]))^(1/6), 0., 0.,
0.}, {0.27712812921102037/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6), 0.,
0., -((2.309401076758503 (-2.625 + \[Delta]))/(-(3 +
2 \[Delta])^2 (-3 + 4 \[Delta]))^(1/6)), 0., 0.}, {0.,
0.27712812921102037/(-(3 + 2 \[Delta])^2 (-3 + 4 \[Delta]))^(
1/6), 0., 0., (
5.923613761885561 -
0.5773502691896257 \[Delta])/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6), 0.}, {0., 0., 0., 0., 0., (
5.7850496972800505 +
1.1547005383792515 \[Delta])/(-(3 + 2 \[Delta])^2 (-3 +
4 \[Delta]))^(1/6)}};


If I plot its eigenvalues using the following code

p = Plot[Evaluate[Eigenvalues[d[\[Delta]]]], {\[Delta], -0.3, 0.3}]


I get the right results as

But I'd like to plot the eigenvalues using another method, constructing the appropriate transformation that afterward, I have the eigenvalues at the diagonal elements of the matrix. If I use the eigenvectors command, I will find such a matrix in which each row is an eigenvector. Then I try to extract data using the following code

For[\[Delta] = -0.3, \[Delta] <= 0.3, \[Delta] = \[Delta] + 0.01,

A = Evaluate[Eigenvectors[d[\[Delta]]]];
MT = A.d[\[Delta]].Inverse[A];

AppendTo[data[[1]], {\[Delta], MT[[1, 1]]}];
AppendTo[data[[2]], {\[Delta], MT[[2, 2]]}];
AppendTo[data[[3]], {\[Delta], MT[[3, 3]]}];
AppendTo[data[[4]], {\[Delta], MT[[4, 4]]}];
AppendTo[data[[5]], {\[Delta], MT[[5, 5]]}];
AppendTo[data[[6]], {\[Delta], MT[[6, 6]]}];
];

ListLinePlot[
{data[[1]], data[[2]], data[[3]], data[[4]], data[[5]], data[[6]]}
]


but the final results are different and wrong:

One can see that after each crossing of the lines, the color of lines changes in a wring way. How can I use the second method in the right way?

As was noted by @HenrikSchumacher the first methods works in the desired way because the eigenvalues of the matrix d are evaluated "analytically", i.e. each eigenvalue is given as a function of the parameter \[Delta]. Therefore, the plot shows them as smooth functions of \[Delta]. In the second plot, where the eigenvalues are calculated from the eigenvectors numerically, the colouring appears to be "wrong" because Mathematica at each parameter value reorders the eigenvectors according to the magnitude of the corresponding eigenvalues, and the smooth dependence on the parameter \[Delta] is lost. In many applications this is considered to be "right", but there is no "right" or "wrong", it depends on what you want.

The second approach could be made to work in the way desired, if Mathematica was able to also evaluate the eigenvectors analytically, i.e. as a function of the parameter \[Delta]. However, unfortunately, the build in algorithm of Mathematica does not do that, it rather reports the error "Not all eigenvectors can be calculated". So for reasons unknown to me the analytical determination of the eigenvectors appears to be more "difficult" for Mathematica than the analytical determination of the eigenvalues.

Here, I provide an implementation of an algorithm for the eigenvectors which works as desired, i.e. calculates the eigenvectors "analytically". I will comment on the algorithm further below. The algorithm relies on the Mathematica implementation for the eigenvalues.

 dim = Length[d[x]];
\[Lambda]M = Eigenvalues /@
Table[Delete[Transpose[Delete[d[x], j]], j], {j, dim}];
\[Lambda] = Eigenvalues[d[x]];
left = Table[\[Lambda][[i]] - \[Lambda][[k]], {i, dim}, {k, dim}] +IdentityMatrix[dim];
left = Times @@ # & /@ left;
right = Table[Product[\[Lambda][[i]] - \[Lambda]M[[j, k]], {k, dim - 1}], {i,dim}, {j, dim}];
v = Table[right[[i, j]]/left[[i]], {i, dim}, {j, dim}];


v contains the norm squared of the elements of the eigenvectors. So, unfortunately, this primitive implementation of the said algorithm does not provide the phases. One can determine the phases as well, but for lack of time I refrain from doing this here, and instead use the numerically calculated eigenvectors of Mathematica to fix the phases. This obviously is a crutch, but sufficient for the demonstration here.

Now one can obtain the desired plot using a similar code like the one used for the second plot, i.e. the eigenvalues determined from the eigenvectors. But now the colouring of the lines is as desired.

len = {}; zeros = Table[0, dim];
data = Table[{\[Delta], old = Eigenvectors[d[\[Delta]]];
new = Sqrt[v /. x -> \[Delta] // Chop];
ord = Flatten[Table[
Position[Table[Abs[old[[i]]] - new[[j]] // Chop, {i, dim}], zeros], {j, dim}]];
len = Join[len, {Length[ord]}];
A = old[[ord]];
Diagonal[A.d[\[Delta]].Transpose[A]]}, {\[Delta], -0.3, .3, .01}];
data = data[[Flatten[Position[len, dim]]]];
ListLinePlot[Table[{#[[1]], #[[2, i]]} & /@ data, {i, Length[A]}]]


The eigenvectors used for the transformation of the matrix d are called A as in the OPs code. Most of the code above serves to implement the mentioned crutch.

Finally, a word about the algorithm used to determine the eigenvectors: This algorithm was rather recently proposed in the paper "Eigenvectors from eigenvalues" https://arxiv.org/pdf/1908.03795.pdf which is coauthored by the fields medal winner Terence Tao. According to this paper the method only works for Hermitian matrices, which is the case for the matrix under consideration here.

• This is awesome! If I can clarify, with this process you are computing the Eigenvectors from the Eigenvalues, as in Tao’s work, then using the numerically computed Eigenvectors from Mathematica to determine the phase shift of the “Eigenvectors from the Eigenvalues”? – CA Trevillian Nov 29 '19 at 19:24
• Not phase shift, here phase is just a fancy word for sign. The algorithm as programmed cannot determine the sign of the components of the eigenvectors. And my crutch fixes the signs. In the cited paper there is also mentioned a facility to determine the sign which I did not yet implement. For more general hermitian matrices than the simple example discussed here, one needs to determine a phase factor Exp[[phi]]. – Michael Weyrauch Nov 29 '19 at 21:05
• Aha! Valuable clarification, thank you! So the signs will be fixed and we will leave with the ordering in-tact? – CA Trevillian Nov 29 '19 at 21:10
• Well, the absolute values of the components are determined from new=Sqrt[...], but I do not know the correct signs. This is then "fixed" by comparison with the numerically determind eigenvector. Clearly "quick and dirty", but I did not have time to properly implement the phase issues as discussed in the cited paper. In a proper implementation one would not need this comparison. – Michael Weyrauch Nov 29 '19 at 21:20
• "here phase is just a fancy word for sign" - since Hermitian matrices are being considered, "phase" is certainly a more general word than "sign", since there is an indeterminacy of the $\exp(i\phi)$ factor for the eigenvector component derived through that procedure. (Additionally, Parlett gives a variant where one does not need to compute eigenvalues of submatrices.) – J. M.'s ennui Nov 30 '19 at 6:07

The issue is the Evaluate in

Plot[Evaluate[Eigenvalues[d[\[Delta]]]], {\[Delta], -0.3, 0.3}]


You certainly use this trick to get the coloring in Plot correct. However, with Evaluate[Eigenvalues[d[\[Delta]]]], you enforce symbolic computation of eigenvalues. The correct ordering for symbolic eigenvalues cannot be figured out until numerical values for \[Delta] are inserted. This is why the ordering of the symbolic eigenvalues may differ from the one of numerically computed eigenvalues.

So a simple fix (and the only that I can provide at the moment) could be to generate a matrix of the plot data and to use ListLinePlot like this (which is basically what you did with the For loop):

data = Transpose@Table[Eigenvalues[d[\[Delta]]], {\[Delta], -0.3, 0.3, 0.01}];
ListLinePlot[data]


# Recoverung continuous eigenvalue curves from the numerical eigenvalues

User CATrevillian asked how we can recover an ordering that makes the eigenvalue functions differentiable. Exploiting that the eigenvalue curves (in this particular case) intersect transversally, one can use the following to obtain this ordering (up to a global permutation).

a = Table[Eigenvalues[d[\[Delta]]], {\[Delta], -0.3, 0.3, 0.01}];
b[[1 ;; 2]] = a[[1 ;; 2]];
Do[
b[[i + 1]] = Flatten[Nearest[a[[i + 1]], b[[i]] + (b[[i]] - b[[i - 1]]), 1]],
{i, 2, Length[a] - 1}]

ListLinePlot[Transpose[b]]


• Dear @Henrik Schumacher, thank you for your comment. Although your solution gives me the same plot, it still has the same problem as my second method and I used. – AYBRXQD Nov 28 '19 at 9:35
• Dear @Henrik Schumacher, even after using Evaluate in computing A, the problem still remains. – AYBRXQD Nov 28 '19 at 10:26
• @AYBRXQD "it still has the same problem as my second method and I used" What I tried to tell you is that the first method is wrong'' in the sense that does not order the eigenvalues by magnitude. – Henrik Schumacher Nov 28 '19 at 11:31
• @AYBRXQD This issue has been discussed several times on this forum. A good discussion is found here. – Michael Weyrauch Nov 28 '19 at 17:24
• @CATrevillian See the recent edit of my post for a quick and dirty approach to recover an ordering that makes the eigenvalue functions differentiable. – Henrik Schumacher Nov 30 '19 at 8:02