对于阶数大于18的情况，Numpy矩阵求逆会给出错误的结果

Question

我在用NumPy反转矩阵时遇到一个问题。奇怪的是，它只给出了直到18阶的正确结果。一旦阶数大于18，就会给出错误的结果。

import numpy as np
from decimal import Decimal
import numpy.matlib
I_1=np.matlib.eye(ngrid,ngrid,k=0,dtype=Decimal)
I_2=np.matlib.eye(ngrid,ngrid,k=1,dtype=Decimal)
I_3=np.matlib.eye(ngrid,ngrid,k=2,dtype=Decimal)
B=I_1 + 10.*I_2 + I_3
B=np.divide(B,12.)
B_inv=np.linalg.inv(B)
print B_inv
C=B.dot(B_inv)
print C

最后一行是为了检查它是否给出了正确的结果。

Answer 1

从技术上讲，您的代码没有任何问题。然而，在数值分析中，你会遇到一个非常基本的问题。这里有两种效果：

1. 浮点错误
2. 问题的条件号

我不会在这里解释它们，因为几乎任何关于数值分析的入门书(例如Kendall Atkinson(都比我在这里做得更好。我只给你留下这个代码片段：

import numpy as np
NGRID = 8

def HilbertMatrix(N, dtype=np.float64):
arr = np.zeros((N,N), dtype=dtype)
for i in range(N):
for j in range(N):
arr[i,j] = 1 / (i + j + 1)
return arr

H = HilbertMatrix(NGRID)
J = np.linalg.inv(H)
C = np.dot(H, J)
print(np.allclose(C, np.eye(NGRID)))
print(np.linalg.norm(C))

试着玩NGRID，看看C离身份矩阵有多远。

Answer 2

这是一个数值分析问题。

由于某种原因，您在numpy数组中保留了十进制对象，这将被转换为Object dtype。当你试图执行矩阵求逆时，矩阵会被转换为浮点(我猜精度很低(。

矩阵的条件数有助于我们理解反转矩阵的难度(有关更正式的定义，请参阅维基百科(。如果您使用np.linal.cond检查B(确保首先转换为numpy-dtype(，您将看到条件号大约为1e18(这是非常非常高的(。

如果你知道一些线性代数，你可以用另一种方法来思考它——矩阵的行列式可以让我们确定矩阵是否奇异。Det(B(=0表示矩阵是奇异的。在你的例子中，你创建了一个矩阵，计算行列式非常容易(上三角(。

矩阵是上三角的，行列式只是对角线的乘积，对角线是(1/12(**(n_grid(。这就是为什么随着ngrid的增长，矩阵几乎是奇异的，对于浮点精度(我猜在引擎盖下会转换为16(，矩阵本质上是奇异/不可逆的。

如果你修复了你的代码，我冒昧地做了这件事，那么隐含的np-dtype是float64，反转可以根据需要工作。

我随意修改了您的代码以使其正常工作(没有定义ngrid((将ngrid更改为19以进行检查(。

import numpy as np
ngrid = 19
I_1 = np.eye(ngrid, ngrid, k=0)
I_2 = np.eye(ngrid, ngrid, k=1)
I_3 = np.eye(ngrid, ngrid, k=2)
B = I_1 + 10.*I_2 + I_3
B = np.divide(B, 12.)
print(np.linalg.cond(B))
B_inv = np.linalg.inv(B)
C = B.dot(B_inv)
assert np.diag(C).sum() == ngrid