python - 使用 scipy.sparse.bmat 从子 block 创建非常大的稀疏矩阵时出错

Question

我正在尝试使用 scipy.sparse.bmat 创建一个矩阵从较小的 csr矩阵 - 我的函数调用在这里:sparse.bmat(HList, format='csr') .生成的矩阵将是一个具有约 26 亿列/行的方阵。但是，当我尝试构造此矩阵时出现以下错误:

    Traceback (most recent call last):
      [...]/lib/python3.7/site-packages/scipy/sparse/construct.py", line 623, in bmat
        return coo_matrix((data, (row, col)), shape=shape).asformat(format)
      [...]/lib/python3.7/site-packages/scipy/sparse/coo.py", line 192, in __init__
        self._check()
      [...]/lib/python3.7/site-packages/scipy/sparse/coo.py", line 283, in _check
        raise ValueError('negative row index found')
    ValueError: negative row index found 

组合矩阵转换为 coo时出现问题格式。我相信这个问题与索引溢出有关，因为完整矩阵的索引不适合 32 位格式(26 亿 > 2^31)。我已经针对这个问题的较小版本测试了我的矩阵构造脚本，并且它工作正常。
This post有一个与我非常相似的问题，但是，那里列出的解决方案对我的情况不起作用。运行那里描述的测试，

>>> scipy.sparse.sputils.get_index_dtype((np.array(10**10),))
<class 'numpy.int64'>

我可以确认 numpy 正在使用 64 位索引。
我的程序是否有其他部分导致溢出？或者问题的根源完全是别的什么？
任何帮助是极大的赞赏!

Answer 1

import numpy as np
from scipy.sparse import coo_matrix, csr_matrix, bmat

a = coo_matrix(([1], ([int(1e9)], [int(1e9)])))
blocks = [a.copy() for i in range(200)]
blocks = [blocks for i in range(20)]

arr = bmat(blocks, format='coo')

首先是第一件事 - 这绝对是可重复的(我使用的是 COO 数组，因为我不想分配 1e11 indptr 数组)。

ValueError: negative row index found

它也无助于转换 a数组索引从 int32 到 int64。实际上，看起来问题完全出在 bmat function 内部。

# convert everything to COO format
for i in range(M):
    for j in range(N):
        if blocks[i,j] is not None:
            A = coo_matrix(blocks[i,j])

首先，它将所有块转换为 COO 矩阵。如果行和列索引适合 int32s，它将使用 int32s(我假设您的索引可以)。稍后它通过添加偏移量(基于块的位置)来计算新的行值。不幸的是，这是它溢出的地方:

for i, j in zip(ii, jj):
    B = blocks[i, j]
    ...
    row[idx] = B.row + row_offsets[i]
    col[idx] = B.col + col_offsets[j]

>>> blocks[2, 0].row
array([1000000000], dtype=int32)

>>> blocks[2, 0].row + 2000000002
array([-1294967294], dtype=int32)

由于溢出(并且因为它在 bmat 中的代码中，您无法从外部访问)，这是一个 scipy 错误。也就是说，如果您复制 scipy bmat 函数并重新键入块索引数组，则可以非常简单地修复它，如下所示:

for i, j in zip(ii, jj):
    B = blocks[i, j]
    ...
    row[idx] = B.row.astype(idx_dtype) + row_offsets[i]
    col[idx] = B.col.astype(idx_dtype) + col_offsets[j]