如何从 scipy 稀疏块矩阵中取回块？

Posted 2023-03-12

【中文标题】如何从 scipy 稀疏块矩阵中取回块？

【英文标题】：How to get the blocks back from a scipy sparse block matrix?

【发布时间】：2019-05-03 14:00:42

【问题描述】：

经过一些矢量化计算，我得到了一个稀疏块矩阵，所有结果都堆叠在相同大小的块中。

>>> A = [[1, 1],
...      [1, 1]]
>>> B = [[2, 2],
...      [2, 2]]
>>> C = [[3, 3],
...      [3, 3]]
>>> results = scipy.sparse.block_diag(A, B, C)
>>> print(results.toarray())
[[1 1 0 0 0 0]
 [1 1 0 0 0 0]
 [0 0 2 2 0 0]
 [0 0 2 2 0 0]
 [0 0 0 0 3 3]
 [0 0 0 0 3 3]]

如果需要，我如何通过提供它们的形状 (2,2) 以有效的方式取回这些数组 A、B、C？

【问题讨论】：

Related.

谢谢！这是一种可行的方法，但它没有利用稀疏性。

【参考方案1】：

In [177]: >>> A = [[1, 1],
     ...: ...      [1, 1]]
     ...: >>> B = [[2, 2],
     ...: ...      [2, 2]]
     ...: >>> C = [[3, 3],
     ...: ...      [3, 3]]
     ...: >>> results = sparse.block_diag([A, B, C])
     ...:      
In [178]: results
Out[178]: 
<6x6 sparse matrix of type '<class 'numpy.int64'>'
    with 12 stored elements in COOrdinate format>

block_diag 不保留输入；而是创建coo 格式矩阵，代表整个矩阵，而不是碎片。

In [194]: results.data
Out[194]: array([1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3], dtype=int64)
In [195]: results.row
Out[195]: array([0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5], dtype=int32)
In [196]: results.col
Out[196]: array([0, 1, 0, 1, 2, 3, 2, 3, 4, 5, 4, 5], dtype=int32)


In [179]: results.A
Out[179]: 
array([[1, 1, 0, 0, 0, 0],
       [1, 1, 0, 0, 0, 0],
       [0, 0, 2, 2, 0, 0],
       [0, 0, 2, 2, 0, 0],
       [0, 0, 0, 0, 3, 3],
       [0, 0, 0, 0, 3, 3]], dtype=int64)

block_diag 将数组传递给sparse.bmat。这反过来从每个矩阵中生成一个coo 矩阵，然后将coo 属性合并为 3 个数组，这些数组是全局稀疏矩阵的输入。

---

还有另一种稀疏格式 bsr 可以保留块（直到转换为 csr 进行计算），但我必须进行实验才能看到这种情况。

让我们从results coo 中创建一个bsr：

In [186]: bresults = sparse.bsr_matrix(results)
In [187]: bresults
Out[187]: 
<6x6 sparse matrix of type '<class 'numpy.int64'>'
    with 12 stored elements (blocksize = 2x2) in Block Sparse Row format>
In [188]: bresults.blocksize
Out[188]: (2, 2)
In [189]: bresults.data
Out[189]: 
array([[[1, 1],
        [1, 1]],

       [[2, 2],
        [2, 2]],

       [[3, 3],
        [3, 3]]], dtype=int64)

所以它推断出有块，就像你想要的那样。

In [191]: bresults.indices
Out[191]: array([0, 1, 2], dtype=int32)
In [192]: bresults.indptr
Out[192]: array([0, 1, 2, 3], dtype=int32)

所以这是一个类似于csr 的存储，但data 是按块排列的。

也许可以在没有block_diag 中介的情况下从您的A,B,C 构造它，但我必须更多地查看文档。

【讨论】：

很好的答案！就我而言，它使我能够以 csr 格式进行所有计算，然后将其转换回 bsr 格式，从那里很好地推断出块，我可以在 bresults.data 中将其取回。非常感谢，我不知道转换为 bsr 格式会隐式恢复结构。

【参考方案2】：

这是一个有趣的小问题。

我不认为有一个函数可以在一行中解决这个问题，但有一种方法可以通过编程来解决。

查看 res.data 打印出来的内容，我在这里使用它。

这适用于所有形状相同的情况。

from scipy.sparse import block_diag

a = [[1, 2, 4],
    [3, 4, 4]]
b = [[2, 2, 1],
    [2, 2, 1]]
c = [[3, 3, 6],
    [3, 3, 6]]

res = block_diag((a, b, c))

def goBack(res, shape):
    s = shape[0]*shape[1]
    num = int(len(res.data)/s)
    for i in range (num):
        mat = res.data[i*s:(i+1)*s].reshape(shape)
        print(mat)

goBack(res, [2,3])

输出：

[[1 2 4]
 [3 4 4]]
[[2 2 1]
 [2 2 1]]
[[3 3 6]
 [3 3 6]]

编辑：

好的，当提供的矩阵的任何元素为零时，这不起作用，因为它不会被计入 res.data。

另外，算了，cleb提供的链接应该对你有帮助。

【讨论】：

这假设了一个特定的稀疏矩阵表示而不指定假设，它假设块是 100% 密集的，没有零。