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

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）都比我在这里做得更好。我只会给你留下这段代码 sn-p：

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.linalg.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