scipy.linalg.block_diag是否有反向操作

Question

我知道scipy提供了一个方便的函数，可以将数组列表转换为块对角矩阵，例如

>>> from scipy.linalg import block_diag
>>> A = [[1, 0],
[0, 1]]
>>> B = [[3, 4, 5],
[6, 7, 8]]
>>> C = [[7]]
>>> D = block_diag(A, B, C)
>>> D
array([[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 3, 4, 5, 0],
[0, 0, 6, 7, 8, 0],
[0, 0, 0, 0, 0, 7]])

这有反向操作吗？也就是说，取一个块对角矩阵和一个块大小列表，然后分解为一个数组列表。

a, b, c = foo(D, block_sizes=[(2,2), (2,3), (1,1)])

如果没有方便的内置方法来实现这一点，那么有没有比在输入数组上循环更好(更高性能(的实现呢？天真的实现可能看起来像这样：

def foo(matrix, block_sizes):
result = []
curr_row, curr_col = 0, 0 
for nrows, ncols in block_sizes:
result.append(matrix[curr_row:curr_row + nrows, curr_col:curr_col + ncols])
curr_row += nrows
curr_col += ncols
return result

Answer 1

查看block_diag:的代码

核心部分分配一个out数组，然后迭代地将块分配给它：

out = np.zeros(np.sum(shapes, axis=0), dtype=out_dtype)
r, c = 0, 0
for i, (rr, cc) in enumerate(shapes):
out[r:r + rr, c:c + cc] = arrs[i]
r += rr
c += cc

你在想象更多的东西吗；表演性的"；？

您想要的是一个形状不同的数组列表。如果没有某种(python级别(的迭代，你将如何获得这些？

查看np.split(及其变体(的代码。它也进行迭代，采用乘法切片来生成数组列表。

Answer 2

FWIW，这里有一些代码使用NumPy来获得非零块索引，即block_diag创建的块所在的相应索引。

import numpy as np
from scipy.linalg import block_diag
def multidim_cumsum_tlbr(a):
# thanks to https://stackoverflow.com/questions/50463366
out = a.cumsum(-1)
for i in range(2, a.ndim+1):
np.cumsum(out, axis=-i, out=out)
return out
def get_block_indices(m):
"""Given a 2D square NumPy array `m` of length `n`, calculates the non-zero
block indices, such that each block contains some non-zero value. All-zero
blocks will yield multiple 1-size blocks.
:param m: a 2D NumPy array
:return: an array of size `[2, n]`, the first row specifying block start
indices, the second row being the end indices (non-inclusive).
>>> get_block_indices(np.array([
>>>     [1,0,0,0,0,0,0,0,0,0,],
>>>     [0,1,0,0,0,0,0,0,0,0,],
>>>     [0,1,1,0,0,0,0,0,0,0,],
>>>     [0,0,0,0,0,0,0,0,0,0,],
>>>     [0,0,0,0,0,0,0,0,0,0,],
>>>     [0,0,0,0,0,0,0,0,0,1,],
>>>     [0,0,0,0,0,0,0,0,0,0,],
>>>     [0,0,0,0,0,0,0,1,0,0,],
>>>     [0,0,0,0,0,0,1,0,0,0,],
>>>     [0,0,0,0,0,0,0,0,0,0,],
>>> ]))
array([[ 0,  1,  3,  4,  5],
[ 1,  3,  4,  5, 10]], dtype=int64)
"""
s = np.shape(m)
if len(s) != 2 or s[0] != s[1]:
raise ValueError('Only accepts 2D NumPy square arrays')

# Indicate non-zero values.
x = (m != 0).astype(int)
# Transpose makes sure that when we count non-zero values for the block,
# we do so for symmetric matrix.
x += x.T
# Count non-zero values along each axis.
x = np.fliplr(multidim_cumsum_tlbr(np.fliplr(x)))
x = np.triu(x)
# Indicate where cumsum is zero.
x = x > 0
# Flag elements where all subsequent elements, in both
# axes, are zero. These are the block ends.
x = ~x[1:,:-1] & ~x[:-1,1:]
# Find block indices.
x = np.nonzero(np.diag(x, 0))[0]
# Calculate end indices and add the final block as a block index.
x = np.hstack([x + 1, [len(m)]])
# Calculate start indices and return the pairs.
x = np.vstack([np.hstack([[0], x[:-1]]), x])
return x

blocks = [
[
[1,0,],
[1,1,],
],
[
[0,0,1,],
[0,0,0,],
[0,0,0,],
],
[
[0,],
],
[
[0,],
],
[
[0,0,0,],
[0,0,0,],
[1,0,0,],
],
[
[0,1,],
[0,0,],
],
[
[-1,1,0,],
[0,-1,1,],
[1,-1,0,],
]
]
m = block_diag(*blocks)
x = get_block_indices(m)
for i, (s, e) in enumerate(x.T):
assert np.all(m[s:e, s:e] == np.array(blocks[i])), (m[s:e, s:e], np.array(blocks[i]))

如果有更好的方式写这篇文章，我不会感到惊讶！我希望人们能改进我的回答。也许还可以处理所有的零块。