【发布时间】:2017-06-23 18:04:59
【问题描述】:
假设我定义了一个大的二次矩阵(例如 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是检查浮点相等性的便捷工具。 -
实际上,一个可能的问题是您的矩阵在构造上几乎是奇异的。我想,反转一个几乎奇异的矩阵会使算法变得非常困难。通常,在这种情况下,机器算术的不可避免的小错误比“正常”操作数更容易累积和爆发。你为什么不尝试一个更乖的例子呢?
标签: python numpy matrix matrix-inverse