假设我定义了一个大的二次矩阵(例如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 中计算出稀疏矩阵的正确逆矩阵?
np.allclose
是用于检查浮点数相等性的方便工具。 - hpaulj