Scipy 简单矩阵乘法的稀疏内存爆炸

Question

我注意到 Scipy 在进行矩阵乘法时必须存储一些中间数组。我认为这在某些情况下会有所帮助，但有时会很痛苦。考虑以下示例：

from scipy.sparse import coo_matrix
n = 100000000000
row  = np.array([0, 0])
col  = np.array([0, n-1])
data = np.array([1, 1])
A = coo_matrix((data, (row, col)), shape=(2, n))

是的，这是一个非常大的矩阵。但是它只有两个非零值。 B = A.dot(A.T) 的结果可以手动求值，因为它只有一个非零值。矩阵 B 使得 B[0, 0] = 2 并且其他地方为零。当我尝试运行时出现内存错误，您可以在下面看到回溯：

    ---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
<ipython-input-32-3e0d3e3c3f13> in <module>
----> 1 A.dot(A.T)

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/base.py in dot(self, other)
    361 
    362         """
--> 363         return self * other
    364 
    365     def power(self, n, dtype=None):

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/base.py in __mul__(self, other)
    479             if self.shape[1] != other.shape[0]:
    480                 raise ValueError('dimension mismatch')
--> 481             return self._mul_sparse_matrix(other)
    482 
    483         # If it's a list or whatever, treat it like a matrix

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/base.py in _mul_sparse_matrix(self, other)
    538 
    539     def _mul_sparse_matrix(self, other):
--> 540         return self.tocsr()._mul_sparse_matrix(other)
    541 
    542     def __rmul__(self, other):  # other * self

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/compressed.py in _mul_sparse_matrix(self, other)
    494 
    495         major_axis = self._swap((M, N))[0]
--> 496         other = self.__class__(other)  # convert to this format
    497 
    498         idx_dtype = get_index_dtype((self.indptr, self.indices,

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/compressed.py in __init__(self, arg1, shape, dtype, copy)
     35                 arg1 = arg1.copy()
     36             else:
---> 37                 arg1 = arg1.asformat(self.format)
     38             self._set_self(arg1)
     39 

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/base.py in asformat(self, format, copy)
    324             # Forward the copy kwarg, if it's accepted.
    325             try:
--> 326                 return convert_method(copy=copy)
    327             except TypeError:
    328                 return convert_method()

~/anaconda3/envs/tfx/lib/python3.7/site-packages/scipy/sparse/coo.py in tocsr(self, copy)
    402             col = self.col.astype(idx_dtype, copy=False)
    403 
--> 404             indptr = np.empty(M + 1, dtype=idx_dtype)
    405             indices = np.empty_like(col, dtype=idx_dtype)
    406             data = np.empty_like(self.data, dtype=upcast(self.dtype))

MemoryError:

输出的是一个2×2的矩阵，所以密不密无所谓。程序试图在如此简单的问题上导致失败的原因是什么？我该如何解决这个问题？

谢谢。

Answer 1

如果你查看源代码，你会发现最终 A.T（A.dot(A.T) 中 dot 的操作数）被转换为 csr_matrix，该矩阵有一个名为 indptr 的属性，其长度为 n+1。使用 n 的值，indptr csr_matrix 将需要数百 GB 的内存。

Answer 2

COO 是一种不适合数学的格式，原因我将留给 reader。

在这种边缘情况下，我认为您最好直接解决问题：

from scipy.sparse import coo_matrix
import numpy as np

n = 100000000000
row  = np.array([0, 0])
col  = np.array([0, n-1])
data = np.array([1, 1])
A = coo_matrix((data, (row, col)), shape=(2, n))

B = A.tocsr()
C = A.tocsr().T

n, m = B.shape[0], C.shape[1]
out_arr = np.zeros((n, m), dtype=A.dtype)

for i in range(n):
    for j in range(m):
        out_arr[i, j] = B[i, :].multiply(C[:, j].T).data.sum()

对于 n 和 m 相当小的任何问题，此解决方法就足够了。

Scipy 简单矩阵乘法的稀疏内存爆炸

Scipy sparse memory explosion with simple matrix multiplication

python

matrix

scipy

sparse-matrix