I was tasked with translating a Mathematica program to Python (Juptyr Notebook) and I'm having trouble getting them to match. The Mathematica program uses the function Eigensystem[] to get the eigenvalues and eigenvectors of a matrix. When I started translating into python, I used the function eig() since it does the same thing. However, I ended up getting different eigenvalues and eigenvectors between the two. The eigenvalues were not off by too much but the eigenvectors ended up being quite different. I had 2 thoughts and was wondering if they are correct:
Just to note, I did notice that the eigenvectors in python are normalized and along the columns while in Mathematica they are typically not and along the rows so that is not the reason. Here is the matrix:
535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 0
228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0
0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0
0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0
0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0
0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0
0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0
0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0
0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0
0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0
0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0
0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 0
0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0 0 98.7
98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0 0 0
0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0 0 0
0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0 0 0
0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0 0 0
0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0 0 0
0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0 0 0
0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0 0 0
0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0 0 0
0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0 0 0
0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0 0 0
0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8 98.7 0
0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 98.7 535.8 228.0
0 0 0 0 0 0 0 0 0 0 0 0 98.7 0 0 0 0 0 0 0 0 0 0 0 228.0 535.8
The code from python:
EigVal, EigVec = eig(H)
print(EigVal)
Results from python:
[114.72374539 141.10213701 154.80916862 224.05697839 240.17122217
313.58869418 326.37713436 329.32297117 350.05452044 362.92295125
381.1300867 402.40762061 415.23547492 656.36452508 669.19237939
690.4699133 708.67704875 721.54547956 742.27702883 745.22286564
758.01130582 831.42877783 847.54302161 916.79083138 930.49786299
956.87625461]
Code for Mathematica:
ES = Eigensystem[H];
ES[[1]]
Results from Mathematica:
953.936, 932.579, 898.44, 853.757, 801.921, 758.013, 747.697,
739.529, 711.095, 700.995, 677.212, 669.422, 536.572, 535.028,
402.178, 394.388, 370.605, 360.505, 332.071, 323.903, 313.587,
269.679, 217.843, 173.16, 139.021, 117.664
I'm not putting in the eigenvectors since there's a bunch and that would take up too much space.
Cannot reproduce:
The code to get the matrix in Python (3.6.15 and numpy 1.19.4)
v1 = [535.8, 228.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 98.7, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 98.7]
v2 = [228.0, 535.8, 98.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 98.7, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
p=[]
for i in range( 0, 26, 2 ):
p.append( np.roll( v1, i ) )
p.append( np.roll( v2, i ) )
p = np.array( p )
p[0,25]=0
p[25,0]=0
Same for Mathematica (11.1.1)
v1 = {535.8, 228.0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 98.7, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 98.7};
v2 = {228.0, 535.8, 98.7, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 98.7, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
p = {};
For[ i = 0, i <= 24, i += 2,
AppendTo[p, RotateRight[v1, i]];
AppendTo[p, RotateRight[v2, i]];
];
p[[1, 26]] = 0;
p[[26, 1]] = 0;
Eigenvalues from Python
EigVal, EigVec = np.linalg.eig( p )
se = sorted( EigVal )
print( se)
provides
[114.72374538612743, 141.1021370092194, 154.80916861965713, 224.05697838622677, 240.1712221708377, 313.5886941750213, 326.3771343636755, 329.3229711673963, 350.05452044459213, 362.9229512535815, 381.1300866993293, 402.4076206114165, 415.2354749186445, 656.364525081355, 669.1923793885824, 690.4699133006694, 708.6770487464171, 721.5454795554072, 742.2770288326029, 745.2228656363234, 758.0113058249784, 831.4287778291624, 847.5430216137748, 916.7908313803451, 930.4978629907807, 956.8762546138717]
and Mathematica gives
Eigenvalues[p] // Sort
{114.724, 141.102, 154.809, 224.057, 240.171, 313.589, 326.377, 329.323, 350.055, 362.923, 381.13, 402.408, 415.235, 656.365, 669.192, 690.47, 708.677, 721.545, 742.277, 745.223, 758.011, 831.429, 847.543, 916.791, 930.498, 956.876}
....and the relative difference after copying the Mathematica values into WMse
print( ( se - WMse ) / se )
are
[-2.21936506e-06 9.70993227e-07 1.08920976e-06 -9.64655213e-08
9.25051868e-07 -9.75242362e-07 4.11682258e-07 -8.75511464e-08
-1.36994491e-06 -1.34316164e-07 2.27479625e-07 -9.42796717e-07
1.14373331e-06 -7.23559283e-07 5.66935001e-07 -1.25565689e-07
6.87850935e-08 6.64622565e-07 3.88434530e-08 -1.80299992e-07
4.03457015e-07 -2.67215718e-07 2.55016846e-08 -1.83923802e-07
-1.47242917e-07 2.66088609e-07]
negligible.
Hence, the problem cannot be reproduced may have other reasons as e.g.the matrix construction.
Note
alternatively the matrices can be constructed like
a = np.diag( [535.8] * 26 )
b = np.diag( [228.0, 98.7] * 12 + [228.0], 1 )
c = np.diag( [228.0, 98.7] * 12 + [228.0], -1 )
d = np.diag( [98.7] * 13 , 13 )
e = np.diag( [98.7] * 13 , -13 )
p = a + b + c + d + e