为什么 numpy 中的矩阵求逆函数和 scipy returns 大二次矩阵的结果不同？

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 中计算正确的 稀疏矩阵的逆矩阵？

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)

通常情况下，这没有问题，因为错误太小，不易察觉。

但是，对于许多计算而言，有些输入难以处理并且可能会过度拉伸数值算法。这当然适用于矩阵求逆，你很不幸选择了如此困难的输入。

您实际上可以通过查看其奇异值来检查矩阵是否可能 'ill-conditioned'，例如 here。以下是使用您的脚本生成的几个矩阵的矩阵条件数（size=200；well-behaved 矩阵的值更接近 1）

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

切换到 well-behaved 矩阵，您的结果应该会大大改善。

为什么 numpy 中的矩阵求逆函数和 scipy returns 大二次矩阵的结果不同？

Why the matrix inversion function in numpy and scipy returns different results with big quadratic matrices?

python

numpy

matrix

matrix-inverse