【问题标题】:Why the matrix inversion function in numpy and scipy returns different results with big quadratic matrices?为什么 numpy 和 scipy 中的矩阵求逆函数对大二次矩阵返回不同的结果?
【发布时间】: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.invspsolve)。

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))

为了检查 AB 的倒数是否相等,我将对差的平方求和。

# 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


【解决方案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)

通常,这没有问题,因为错误太小而无法注意到。

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

您实际上可以通过查看矩阵的奇异值来检查矩阵是否“病态”,例如here。以下是您的脚本生成的几个矩阵的矩阵条件数(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

切换到表现良好的矩阵,您的结果应该会大大改善。

【讨论】:

    猜你喜欢
    • 2021-02-23
    • 2017-06-07
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-05-11
    • 2021-12-20
    • 1970-01-01
    • 2018-02-17
    相关资源
    最近更新 更多