为什么在numpy和scipy中的矩阵求逆函数在处理大型二次矩阵时返回不同的结果？

Question

假设我定义了一个大的二次矩阵（例如150x150）。一次是numpy数组（矩阵A），一次是scipy稀疏数组（矩阵B）。

import numpy as np
import scipy as sp

from scipy.sparse.linalg import spsolve

size = 150
A = np.zeros((size, size))
length = 1000
# Set some random numbers at random places
A[np.random.randint(0, size, (2, length)).tolist()] = \
    np.random.randint(0, size, (length, ))
B = sp.sparse.csc_matrix(A)

现在我计算两个矩阵的逆矩阵。对于矩阵B，我使用了两种方法来计算逆矩阵(sp.sparse.linalg.inv和spsolve)。

epsilon = 10.**-8 # Is needed to prevent singularity of each matrix

inv_A = np.linalg.pinv(A+np.eye(size)*epsilon)
inv_B = sp.sparse.linalg.inv(B+sp.sparse.identity(size)*epsilon)
inv_B2 = spsolve(B+sp.sparse.identity(size)*epsilon, sp.sparse.identity(size))

为了检查A和B的倒数是否相等，我将对差值进行平方求和。

# Is not equal zero, question: Why?
# Sometimes very small (~+-10**-27), sometimes very big (~+-10**5)
print("np.sum((inv_A - inv_B)**2): {}".format(np.sum((inv_A - inv_B)**2)))
# Should be zero
print("np.sum((inv_B - inv_B2)**2): {}".format(np.sum((inv_B - inv_B2)**2)))

问题在于：如果我使用小矩阵，例如 10x10，那么 numpy 和 scipy 的求逆函数之间的误差非常小（约为 ~+-10**-32）。但是对于大矩阵（例如 500x500），我需要稀疏版本。

我这里有什么做错了吗？还是有可能在 Python 中计算出稀疏矩阵的正确逆矩阵？

- PiMathCLanguage

---

相对误差有多大？ - Paul Panzer

你是指我的问题可能的相对误差还是两个矩阵之间的相对误差？ - PiMathCLanguage

np.allclose是用于检查浮点数相等性的方便工具。 - hpaulj

实际上，一个可能的问题是你所构建的矩阵接近奇异。我猜，对于一个几乎奇异的矩阵来说，求逆操作会让算法变得非常困难。通常，在这种情况下，机器计算中不可避免的小误差会比“正常”的运算更容易积累并且放大。为什么不尝试一个更规范的例子呢？ - Paul Panzer

您IP地址为143.198.54.68，由于运营成本限制，当前对于免费用户的使用频率限制为每个IP每72小时10次对话，如需解除限制，请点击左下角设置图标按钮（手机用户先点击左上角菜单按钮）。 - PiMathCLanguage

Answer 1

你的问题的答案是：由于您不幸选择了示例矩阵。让我详细说明一下。

机器精度有限，因此浮点算术很少是100％准确的。只需尝试

>>> np.linspace(0, 0.9, 10)[1:] == np.linspace(0.1, 1, 10)[:-1]
array([ True,  True,  True,  True,  True, False,  True,  True,  True], dtype=bool)

通常情况下，这不是问题，因为错误很小以至于不易察觉。

然而，在许多计算中，有些输入很难处理，可能会过度拉伸数值算法。这当然适用于矩阵求逆，你不幸选择了这样具有挑战性的输入。

你实际上可以通过查看矩阵奇异值来检查一个矩阵是否“病态”，例如在这里看到。以下是使用你的脚本生成的几个矩阵的条件数（size = 200；一个行为良好的矩阵其值更接近于1）

971899214237.0
5.0134186641e+12
36848.0807109
958492416768.0
1.66615247737e+16
1.42435766189e+12
1954.62614384
2.35259324603e+12
5.58292606978e+12

如果您采用良好的矩阵方法，您的结果应该会显著提高。

- Paul Panzer