将大稀疏矩阵转换为COO时出错

Question

我在尝试把两个大的CSR矩阵上下堆叠时遇到了一个问题：

    /usr/lib/python2.7/dist-packages/scipy/sparse/coo.pyc in _check(self)
    229                 raise ValueError('negative row index found')
    230             if self.col.min() < 0:
--> 231                 raise ValueError('negative column index found')
    232
    233     def transpose(self, copy=False):

ValueError: negative column index found

我可以很简单地通过尝试把一个大的lil矩阵转换成coo矩阵来重现这个错误。下面的代码在N=10**9时可以正常运行，但在N=10**10时就失败了。

from scipy import sparse
from numpy import random
N=10**10
x = sparse.lil_matrix( (1,N) )
for _ in xrange(1000):
    x[0,random.randint(0,N-1)]=random.randint(1,100)

y = sparse.coo_matrix(x)

我是不是碰到了coo矩阵的大小限制？有没有什么办法可以解决这个问题？

Answer 1

看起来你遇到了32位整数的限制。这里有个简单的测试：

In [14]: np.array([10**9, 10**10], dtype=np.int64)
Out[14]: array([ 1000000000, 10000000000])

In [15]: np.array([10**9, 10**10], dtype=np.int32)
Out[15]: array([1000000000, 1410065408], dtype=int32)

目前，大多数稀疏矩阵的表示方法都假设使用32位整数作为索引，所以它们无法支持那么大的矩阵。

补充说明：从0.14版本开始，scipy现在支持64位索引。如果你能升级一下，这个问题就能解决了。

Answer 2

有趣的是，你的第二个例子在我的安装环境中运行得很好。

错误信息“发现负列索引”听起来像是某个地方出现了溢出。我检查了最新的源代码，得到了以下结果：

• 实际的索引数据类型是在 scipy.sparse.sputils.get_index_dtype 中计算的。
• 这个错误信息来自模块 scipy.sparse.coo。

这个异常是由以下类型的代码引起的：

    idx_dtype = get_index_dtype(maxval=max(self.shape))
    self.row = np.asarray(self.row, dtype=idx_dtype)
    self.col = np.asarray(self.col, dtype=idx_dtype)
    self.data = to_native(self.data)

    if nnz > 0:
        if self.row.max() >= self.shape[0]:
            raise ValueError('row index exceeds matrix dimensions')
        if self.col.max() >= self.shape[1]:
            raise ValueError('column index exceeds matrix dimensions')
        if self.row.min() < 0:
            raise ValueError('negative row index found')
        if self.col.min() < 0:
            raise ValueError('negative column index found')

这明显是一个溢出错误，可能是在 - 2 的 31 次方。

如果你想调试这个问题，可以试试：

import scipy.sparse.sputils
import numpy as np

scipy.sparse.sputils.get_index_dtype((np.array(10**10),))

它应该返回 int64。如果没有返回，那问题就出在这里。

你使用的是哪个版本的 SciPy？