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

Question

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

解释了如何对单个矩阵执行此操作，但我有一大堆不同的矩阵需要制作成块矩阵。

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

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

Answer 1

对于未来的读者：我找到了一种完全在 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 对应函数慢得多，可能是因为需要对索引的 concatenation/shift 使用扫描。欢迎任何关于如何解决这个问题的想法！

Answer 2

谢谢Tool

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

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

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

Build block diagonal matrix from many small matrices in numpy / theano

python

arrays

numpy

theano