np.linalg.inv使用疑问：求逆后点乘原矩阵出现浮点精度误差

代码中定义矩阵数组a及其理论逆矩阵b，二者点乘应得到单位矩阵I：手动定义的b与a点乘确实输出标准单位矩阵，结果符合预期；但调用np.linalg.inv计算a的逆矩阵得到b1时，b1数值和b几乎一致，但a与b1点乘得到的I1和标准单位矩阵存在极小偏差。已确认矩阵a行列式非零，不属于奇异矩阵。
示例代码如下：

import numpy as np

a = np.array([[0,1,2], [1,2,3], [3,1,1]])
b = np.array([[1/2,-1/2,1/2], [-4,3,-1],[5/2,-3/2,1/2]])

I = np.dot(a,b)
print(I)
# 输出
# [[1., 0., 0.],
#  [0., 1., 0.],
#  [0., 0., 1.]]

b1 = np.linalg.inv(a)
I1 = np.dot(a,b1)
print(I1)
# 输出
# [[1.00000000e+00, 0.00000000e+00, 1.11022302e-16],
#  [8.88178420e-16, 1.00000000e+00, 2.22044605e-16],
#  [0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]

该偏差是双精度浮点数的固有精度限制导致的，不属于函数计算错误：

• 手动定义的b是精确分数值转换为64位浮点数，和a做矩阵乘法时，每一步乘加运算的结果刚好能被浮点数精确表示，因此输出完全标准的单位矩阵，没有可见误差。
• np.linalg.inv底层通过LU分解等数值线性代数算法迭代求解逆矩阵，计算过程中会产生多次浮点舍入。观测到的非零偏差量级在1e-16，刚好是双精度浮点数的机器精度（即能表示的最小相对误差）级别，属于数值计算中的正常现象。

• 绝大多数常规数值计算场景下，1e-16量级的误差不会对后续结果产生实质影响，不需要额外处理。
• 若需要判断矩阵乘积是否为单位矩阵，禁止使用==做精确值比较，应使用numpy自带的带公差比较函数：例如调用np.allclose(I1, np.eye(a.shape[0]))即可在忽略浮点误差的前提下判断结果是否符合预期，该函数默认的相对误差、绝对误差阈值完全覆盖常规浮点计算的误差范围。
• 若有极高精度计算需求，可以将数组的数据类型指定为更高精度的浮点数（如np.float128，该类型的支持依赖运行环境的硬件与操作系统），但会同步提升内存占用和计算耗时，普通场景无必要使用。

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