为啥 scipy sparse 和 numpy 数组的乘法函数会给出不同的结果？

Posted 2023-02-25

【中文标题】为啥 scipy sparse 和 numpy 数组的乘法函数会给出不同的结果？

【英文标题】：Why multiplication functions of scipy sparse and numpy arrays give different results?为什么 scipy sparse 和 numpy 数组的乘法函数会给出不同的结果？

【发布时间】：2020-01-13 09:10:57

【问题描述】：

我在 Python 2.7 中有两个矩阵：一个密集 A_dense 和另一个稀疏矩阵 A_sparse。我对计算元素乘法和总和感兴趣。有两种方法可以做到这一点：使用numpy's multiplication 或scipy sparse multiplication。我希望他们给出完全相同的结果，但执行时间不同。但我发现它们对某些矩阵大小给出了不同的结果。

import numpy as np
from scipy import sparse
L=2000
np.random.seed(2)
rand_x=np.random.rand(L)
A_sparse_init=np.diag(rand_x, -1)+np.diag(rand_x, 1)
A_sparse=sparse.csr_matrix(A_sparse_init)
A_dense=np.random.rand(L+1,L+1)
print np.sum(A_sparse.multiply(A_dense))-np.sum(np.multiply(A_dense[A_sparse.nonzero()], A_sparse.data))

输出：

1.1368683772161603e-13

如果我选择 L=2001，那么输出是：

0.0

为了使用两种不同的乘法方法检查差异的大小依赖性，我写道：

L=100
np.random.seed(2)
N_loop=100
multiply_diff_arr=np.zeros(N_loop)
for i in xrange(N_loop):
    rand_x=np.random.rand(L)
    A_sparse_init=np.diag(rand_x, -1)+np.diag(rand_x, 1)
    A_sparse=sparse.csr_matrix(A_sparse_init)
    A_dense=np.random.rand(L+1,L+1)
    multiply_diff_arr[i]=np.sum(A_sparse.multiply(A_dense))-np.sum(np.multiply(A_dense[A_sparse.nonzero()], A_sparse.data))
    L+=1

我得到了以下情节：

谁能帮我理解发生了什么？难道我们不期望两种方法之间的差异至少是 1e-18 而不是 1e-13？

【参考方案1】：

我没有完整的答案，但这可能有助于找到答案：

在后台，scipy.sparse 将转换为 coo 格式并执行以下操作：

ret = self.tocoo()
if self.shape == other.shape:
    data = np.multiply(ret.data, other[ret.row, ret.col])

那么问题是为什么这两个操作会给出不同的结果：

ret = A_sparse.tocoo()
c = np.multiply(ret.data, A_dense[ret.row, ret.col])
ret.data = c.view(type=np.ndarray)

c.sum() - ret.sum()

-1.1368683772161603e-13

编辑：

差异源于不同的默认值在哪个轴上首先add.reduce。 例如：

A_sparse.multiply(A_dense).sum(axis=1).sum()
A_sparse.multiply(A_dense).sum(axis=0).sum()

Numpy 默认为0。

【讨论】：

Why does changing the sum order returns a different result?