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

How to get the blocks back from a scipy sparse block matrix?

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

>>> 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?

这是个有趣的小问题。

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

查看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提供的link应该可以帮到你

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 构建它,但我必须更多地查看文档。