Python/numpy与Mathematica计算矩阵特征值/向量结果不一致排查

我承担了将Mathematica程序移植为Python（Jupyter Notebook）代码的任务，移植过程中发现两侧计算结果无法匹配：Mathematica程序调用Eigensystem[]函数计算矩阵的特征值与特征向量，Python侧选用功能对应的eig()函数实现相同逻辑，但最终计算结果存在明显差异——特征值偏差幅度较小，但特征向量差异较大。
已提前排除特征向量存储与归一化规则的影响：Python返回的特征向量为归一化形式、按列存储，Mathematica返回的特征向量默认不归一化、按行存储，该差异并非导致结果不匹配的原因。

待验证猜想 #

• 待计算矩阵规模为26×26，较大的矩阵尺寸是否会导致两个平台的数值计算结果出现差异？
• Mathematica更擅长符号数学运算、计算过程无舍入误差，这一特性是否是造成结果偏差的原因？

待计算的26×26矩阵 #

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

两侧实现代码与运行结果 #

Python侧 #

实现代码：

EigVal, EigVec = eig(H)
print(EigVal)

运行得到的特征值结果：

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

Mathematica侧 #

实现代码：

ES = Eigensystem[H];
ES[[1]]

运行得到的特征值结果：

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

注：因特征向量数据量较大，此处未附上对应的计算结果。

内容的提问来源于stack exchange，提问作者SlickDickie927