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

Posted 2023-03-11

【中文标题】使用 scipy.sparse.bmat 从子块创建非常大的稀疏矩阵时出错

【英文标题】：Error creating very large sparse matrix from sub-blocks using scipy.sparse.bmat

【发布时间】：2021-10-21 13:53:33

【问题描述】：

我正在尝试使用较小的 csr 矩阵中的 scipy.sparse.bmat 创建一个矩阵 - 我的函数调用在这里：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 位索引。

我的程序是否有其他部分导致溢出？或者问题的根源完全是别的什么？

非常感谢任何帮助！

【问题讨论】：

看起来你走在正确的轨道上，但我不认为我会尝试在我的普通电脑上测试这么大的东西。您可以从源代码中获取bmat 的代码，并通过一些形状和数据类型检查制作您自己的版本，以验证它是否将正确的shape、row、col 传递给coo_matrix。跨度>

get_index_dtype 运行两次，分别在 bmat 和 coo._check 中。

您还可以验证它是否运行：sparse.coo_matrix((int(3e9),int(2.6e9)),int)。以前没见过这样的东西，所以只能根据我在代码中看到的内容提出建议。

@hpaulj 感谢您的建议，原来是问题所在

【参考方案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]

【讨论】：

太棒了，感谢您解决这个问题！我已向 scipy 的开发人员提交了一个描述此错误的问题，因此它得到了修复。