python - 为什么 Scipy 对该矩阵的特征值得到不同的答案？

Question

在 Scipy 中，我有一个 7x7 矩阵，其中一些条目大到 10^22。

[[7.87720699e+21 5.81984000e+12 4.06283195e+18 1.36914426e+15
  6.84296262e+21 1.18920842e+20 2.30532710e+21]
 [5.81984000e+12 6.32869618e+03 3.12159108e+09 8.49240505e+05
  5.05790011e+12 8.25020510e+10 1.71012096e+12]
 [4.06283195e+18 3.12159108e+09 2.10517370e+15 6.90360153e+11
  3.52958137e+18 6.08623401e+16 1.18959479e+18]
 [1.36914426e+15 8.49240505e+05 6.90360153e+11 2.68568446e+08
  1.18907574e+15 2.15022766e+13 3.99729760e+14]
 [6.84296262e+21 5.05790011e+12 3.52958137e+18 1.18907574e+15
  5.94451366e+21 1.03298010e+20 2.00265814e+21]
 [1.18920842e+20 8.25020510e+10 6.08623401e+16 2.15022766e+13
  1.03298010e+20 1.81926760e+18 3.47747090e+19]
 [2.30532710e+21 1.71012096e+12 1.18959479e+18 3.99729760e+14
  2.00265814e+21 3.47747090e+19 6.74706499e+20]]

当我调用spicy.linalg.eigvalsh时，我得到

scipy.linalg.eigvalsh(square_matrix)
array([-1.65239967e+05,  3.55247340e+04,  2.64944833e+06,  2.26542682e+09,\n        2.01743752e+14,  5.56910661e+16,  1.44981926e+22])

这不好，因为矩阵是对称正半定矩阵(其形式为 AA^T)，因此它应该具有所有非负特征值。当我只要求它提供最小特征值时，它得到的答案与上次完全不同。

scipy.linalg.eigvalsh(square_matrix, eigvals = (0,0))
array([-464577.85826165])

这是怎么回事？难道只是数字太大，导致计算出现巨大误差？

感谢您的帮助，我对这样的数值实验有点陌生。

编辑:这是另一个矩阵，以列表形式发布:

[[178484429459288.62, 1.262534539362581e+16, 5756113437609551.0, 1.5274842899247696e+17, 6343247145.960138, 3119899812227.452, 2.6847705451184886e+17], [1.262534539362581e+16, 9.311154944318847e+17, 4.25182007604871e+17, 1.0968899948471126e+19, 418362473621.39294, 236566750583994.28, 1.929426966567391e+19], [5756113437609551.0, 4.25182007604871e+17, 1.9416578714719302e+17, 5.00372260140452e+18, 190262613997.92053, 108142802848362.1, 8.801791711005543e+18], [1.5274842899247696e+17, 1.0968899948471126e+19, 5.00372260140452e+18, 1.3144272477543575e+20, 5288261203200.205, 2737268307603631.5, 2.3109381789197917e+20], [6343247145.960138, 418362473621.39294, 190262613997.92053, 5288261203200.205, 258730.2206841113, 99079200.39025305, 9282686524313.389], [3119899812227.452, 236566750583994.28, 108142802848362.1, 2737268307603631.5, 99079200.39025305, 61292039604.753365, 4817296126969346.0], [2.6847705451184886e+17, 1.929426966567391e+19, 8.801791711005543e+18, 2.3109381789197917e+20, 9282686524313.389, 4817296126969346.0, 4.0629951073761985e+20]]

scipy.linalg.eigvalsh(square_matrix, eigvals = (0,0))
> array([-9362.07065027])

scipy.linalg.eigvalsh(square_matrix)
> array([3.04005337e+01, 1.06920150e+04, 9.94140064e+06, 3.42415026e+09,\n       2.27074109e+13, 1.99821497e+16, 5.38847690e+20])

我已确认我使用的是 Scipy 1.1.0。如果有影响的话，这些是在脚本中间调试期间获得的。使用 VS Code 作为 IDE。

Answer 1

矩阵的条目最多为 1e22，但其最小特征值可能要小得多，可能为 1000 左右。当使用 double 进行复杂计算时(1e-16每一步的相对误差)，尤其是近似特征值的迭代计算，得到大小为1e-12的相对误差并不奇怪。从绝对值来看，它变成了1e10，这使得最小特征值完全相形见绌，并且很容易使其成为负值。

搜索 k 个最小或最大特征值是一项不同的任务:SciPy 不仅仅计算所有特征值并截断输出。在 source根据 eigvals 参数的存在，人们可以看到对基于 LAPACK 的代码的不同调用，并且我确信该参数在 LAPACK 内部也会产生影响。所有这些意味着您得到的随机错误很可能会有所不同，具体取决于参数eigvals。

I'm trying to estimate the condition number

我建议使用np.linalg.cond为此。

一般来说，处理 A 的奇异值比处理 A 的特征值乘以其转置(条件数 gets nearly squared )更好。