从 numpy / theano 中的许多小矩阵构建块对角矩阵答案

【问题标题】：Build block diagonal matrix from many small matrices in numpy / theano从 numpy / theano 中的许多小矩阵构建块对角矩阵
【发布时间】：2018-04-15 20:04:12
【问题描述】：

我想从相同形状的方阵列表中构建块对角矩阵。

This answer 解释了如何对单个矩阵进行处理，但我有一大堆不同的矩阵需要制作成块矩阵。

我想在 GPU 上的 Theano 中执行此操作，因此性能是必需的（并且该功能必须存在于 Theano 中）。

详情： 这样做的原因是当有许多小矩阵（例如，大约 10000 个 7x7 矩阵）时，在 GPU 上加速特征值/向量的计算。我不想分别获取每个小矩阵的特征值（在 GPU 上非常慢），而是想在块对角矩阵上执行大 EVD（与小矩阵相同的特征值）。希望这会更快，并且矩阵的稀疏性不会产生太多开销（或者 EIGH 可能会利用这一点）。

【问题讨论】：

scipy 中有一个块对角方法，也许你可以使用它。为什么需要在theano中构造矩阵？ docs.scipy.org/doc/scipy-1.0.0/reference/generated/…
@JacquesKvam 我注意到 scipy 具有此功能，尽管我正在寻找 numpy / theano 中的实现。为什么要问原因？我正在使用 theano，所以我需要在 theano 中使用它！
如果你只需要构造一次矩阵，你可以用numpy/scipy来做，然后转换成theano。很难说没有更多细节。
@JacquesKvam 添加了一些细节。当然，如果我只需要执行一次，我就不会在 theano 中执行：问题是我必须在训练 DNN 期间执行多次（实际上，作为梯度下降算法的一部分）。
This Mathematica answer does it 虽然我在用 numpy 翻译代码时遇到了麻烦。

标签： python arrays numpy theano

【解决方案1】：

谢谢Tool

我稍微修改了代码，使它在第 k 个对角线上返回 x。

def block_diagonal(x, k):
''' x should be a tensor-3 (#num matrices, n,n)
    k : int
    Diagonal in question. it is 0 in case of main diagonal. 
    Use k>0 for diagonals above the main diagonal, and k<0 for diagonals below the main diagonal.
'''

shape = x.shape
n = shape[-1]

absk = abs(k)

indx = np.repeat(np.arange(n),n)
indy = np.tile(np.arange(n),n)

indx = np.concatenate([indx + a * n for a in range(shape[0])])
indy = np.concatenate([indy + a * n for a in range(shape[0])])

if k<0: 
    indx += n*absk
else:
    indy += n*absk

block = np.zeros(((shape[0]+absk)*n,(shape[0]+absk)*n))
block[(indx,indy)] = x.flatten()

return block

【讨论】：

【解决方案2】：

对于未来的读者：我找到了一种纯粹在 NumPy 中执行此操作的方法，它比 SciPy 函数略快：

def block_diagonal(x):
" x should be a tensor-3 (#num matrices, n,n) "

shape = x.shape
n = shape[-1]

indx = np.repeat(np.arange(n),n)
indy = np.tile(np.arange(n),n)

indx = np.concatenate([indx + k * n for k in range(shape[0])])
indy = np.concatenate([indy + k * n for k in range(shape[0])])

block = np.zeros((shape[0]*n,shape[0]*n))
block[(indx,indy)] = x.flatten()

return block

这个实现只是建立块所在的索引，然后填充它！

时间：

matrix_list = [np.arange(i,i+49).reshape((7,7)) for i in range(1000)]

matlist = np.array(matrix_list)

%timeit block_diagonal(matlist)
25.6 ms ± 145 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit scipy.linalg.block_diag(*matrix_list)
28.6 ms ± 227 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

matrix_list = [np.arange(i,i+49).reshape((7,7)) for i in range(5000)]
matlist = np.array(matrix_list)

%timeit block_diagonal(matlist)
141 ms ± 75.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit scipy.linalg.block_diag(*matrix_list)
157 ms ± 201 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

注意：Theano 中的相同函数比其 NumPy 对应的函数慢得多，可能是因为需要使用扫描来进行索引的连接/移位。欢迎任何关于如何解决这个问题的想法！

【讨论】：