numpy 矩阵求逆的四舍五入误差

Question

我在计算我的BinvA矩阵时，发现(1,1)这个位置的值非常奇怪。
我只是想把B矩阵反转，然后做(B^-1)A的乘法。

我知道手动计算时，(1,1)的位置应该是0，但我得到的却是1.11022302e-16。这个值怎么回事？我知道浮点数不能完全准确地表示，但为什么会给我这么不准确的结果，而不是四舍五入到0呢？有没有什么办法可以让它更准确一些？

这是我的代码：

import numpy as np

A = np.array([[2,2],[4,-1]],np.int)
A = A.transpose()


B = np.array([[1,3],[-1,-1]],np.int)
B = B.transpose()

Binv = np.linalg.inv(B) #calculate the inverse

BinvA = np.dot(Binv,A) 
print(BinvA)

这是我的打印语句：

[[  1.11022302e-16  -2.50000000e+00]
 [ -2.00000000e+00  -6.50000000e+00]]

Answer 1

这不是一个完整的答案，但可能会给你一些启发。你真正想要的是使用小数（Decimal）进行数学运算的numpy数组。你可能会合理地想到尝试：

import numpy as np
from decimal import Decimal
A = np.array([[2,2],[4,-1]],np.int)
for i, a in np.ndenumerate(A):
    A[i] = Decimal(a)
    print type(A[i])

但是，遗憾的是，小数并不是numpy默认支持的数据类型，所以每次你试图把小数放进数组时，它都会把它转换成浮点数（float）。

一种可能的解决办法是这样设置数据类型：

def decimal_array(arr):
    X = np.array(arr, dtype = Decimal)
    for i, x in np.ndenumerate(X): X[i] = Decimal(x)
    return X

A = decimal_array([[2,2],[4,-1]])
B = decimal_array([[1,3],[-1,-1]])

A = A.transpose()
B = B.transpose()
Binv = np.linalg.inv(B) #calculate the inverse

但是现在，如果你

print Binv.dtype

你会发现它又被转换回浮点数了。原因是linalg.inv（和许多其他函数一样）会寻找B的“公共类型”，也就是它认为可以强制转换你数组元素的标量类型。

不过，这也不是完全没有希望。我查了一下，看看是否可以通过创建自定义数据类型来解决这个问题，但结果发现标量（整数、浮点数等）根本不是数据类型。相反，你可能想要做的是注册一个新的标量——也就是小数——正如在关于标量的文章中所说的那样。你会看到一个链接指向Numpy的C-API（别害怕）。在页面上搜索“register”和“scalar”来开始了解。

Answer 2

当你计算逆矩阵时，你的数组会被转换成 float64 类型，这种类型的 机器精度 是 1e-15。这个精度值可以理解为浮点数的 相对量化步长。

如果你对浮点数的数据类型有疑问，可以通过 finfo 函数向 numpy 查询相关信息。在这种情况下：

np.finfo('float64')
finfo(resolution=1e-15, 
      min=-1.7976931348623157e+308, max=1.7976931348623157e+308, 
      dtype=float64)

所以，从技术上讲，你的值如果小于 eps，对于 float64 类型来说，这实际上是非常准确地表示了 0。

如果只是这种表示方式让你觉得困扰，你可以告诉 numpy 不要打印小于等于 0 的浮点数（即 1 个 epsilon 或更小的数），方法是：

np.set_printoptions(suppress=True)

这样一来，你的打印语句返回的结果就是：

[[ 0.  -2.5]
 [-2.  -6.5]]

需要注意的是，这其实是所有浮点数实现中普遍存在的一个数值问题。你可以在 Stack Overflow 上找到更多关于浮点数舍入误差的信息：

• 为什么浮点数不准确？

或者在网上查找：

• 浮点数精度问题
• 每个计算机科学家都应该知道的浮点运算知识