从给定的numpy数组创建块对角线numpy数组
从给定的numpy数组创建块对角线numpy数组
提问于 2015-11-04 04:20:58
我有一个具有相等列数和行数的二维numpy数组。我想把它们排列成一个较大的数组，其中较小的在对角线上。应该可以指定起始矩阵在对角线上的频率。例如：
a = numpy.array([[5, 7],
[6, 3]])所以如果我想让这个数组在对角线上重复两次，那么期望的输出将是：
array([[5, 7, 0, 0],
[6, 3, 0, 0],
[0, 0, 5, 7],
[0, 0, 6, 3]])3次：
array([[5, 7, 0, 0, 0, 0],
[6, 3, 0, 0, 0, 0],
[0, 0, 5, 7, 0, 0],
[0, 0, 6, 3, 0, 0],
[0, 0, 0, 0, 5, 7],
[0, 0, 0, 0, 6, 3]])对于任意大小的起始数组(仍然认为起始数组具有相同的行数和列数)，有没有一种快速的方法来实现numpy方法？
回答 3
Stack Overflow用户
回答已采纳
发布于 2015-11-04 04:23:15
方法#1
numpy.kron的经典案例-
np.kron(np.eye(r,dtype=int),a) # r is number of repeats示例运行-
In [184]: a
Out[184]:
array([[1, 2, 3],
[3, 4, 5]])
In [185]: r = 3 # number of repeats
In [186]: np.kron(np.eye(r,dtype=int),a)
Out[186]:
array([[1, 2, 3, 0, 0, 0, 0, 0, 0],
[3, 4, 5, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 2, 3, 0, 0, 0],
[0, 0, 0, 3, 4, 5, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 2, 3],
[0, 0, 0, 0, 0, 0, 3, 4, 5]])方法#2
另一个使用diagonal-viewed-array-assignment的有效方法--
def repeat_along_diag(a, r):
m,n = a.shape
out = np.zeros((r,m,r,n), dtype=a.dtype)
diag = np.einsum('ijik->ijk',out)
diag[:] = a
return out.reshape(-1,n*r)示例运行-
In [188]: repeat_along_diag(a,3)
Out[188]:
array([[1, 2, 3, 0, 0, 0, 0, 0, 0],
[3, 4, 5, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 2, 3, 0, 0, 0],
[0, 0, 0, 3, 4, 5, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 2, 3],
[0, 0, 0, 0, 0, 0, 3, 4, 5]])Stack Overflow用户
发布于 2018-05-28 20:40:56
import numpy as np
from scipy.linalg import block_diag
a = np.array([[5, 7],
[6, 3]])
n = 3
d = block_diag(*([a] * n))
d
array([[5, 7, 0, 0, 0, 0],
[6, 3, 0, 0, 0, 0],
[0, 0, 5, 7, 0, 0],
[0, 0, 6, 3, 0, 0],
[0, 0, 0, 0, 5, 7],
[0, 0, 0, 0, 6, 3]], dtype=int32)但是看起来对于小的n，np.kron解决方案要快一点。
%timeit np.kron(np.eye(n), a)
12.4 µs ± 95.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit block_diag(*([a] * n))
19.2 µs ± 34.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)然而，对于n= 300，例如，block_diag要快得多：
%timeit block_diag(*([a] * n))
1.11 ms ± 32.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
%timeit np.kron(np.eye(n), a)
4.87 ms ± 31 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)Stack Overflow用户
发布于 2019-09-20 02:51:43
对于矩阵的特殊情况，简单的切片比numpy.kron() (最慢的)快得多，基本上与numpy.einsum()-based方法(来自@Divakar -based)相当。与scipy.linalg.block_diag()相比，对于较小的arr，它的性能更好，这在一定程度上独立于块重复次数。
请注意，block_diag_view()在较小输入上的性能可以很容易地通过Numba进一步提高，但对于较大的输入会得到较差的性能。
使用Numba，完全显式循环和并行化(block_diag_loop_jit())如果重复次数很少，将再次获得与block_diag_einsum()相似的结果。
总体而言，性能最好的解决方案是block_diag_einsum()和block_diag_view()。
import numpy as np
import scipy as sp
import numba as nb
import scipy.linalg
NUM = 4
M = 9
def block_diag_kron(arr, num=NUM):
return np.kron(np.eye(num), arr)
def block_diag_einsum(arr, num=NUM):
rows, cols = arr.shape
result = np.zeros((num, rows, num, cols), dtype=arr.dtype)
diag = np.einsum('ijik->ijk', result)
diag[:] = arr
return result.reshape(rows * num, cols * num)
def block_diag_scipy(arr, num=NUM):
return sp.linalg.block_diag(*([arr] * num))
def block_diag_view(arr, num=NUM):
rows, cols = arr.shape
result = np.zeros((num * rows, num * cols), dtype=arr.dtype)
for k in range(num):
result[k * rows:(k + 1) * rows, k * cols:(k + 1) * cols] = arr
return result
@nb.jit
def block_diag_view_jit(arr, num=NUM):
rows, cols = arr.shape
result = np.zeros((num * rows, num * cols), dtype=arr.dtype)
for k in range(num):
result[k * rows:(k + 1) * rows, k * cols:(k + 1) * cols] = arr
return result
@nb.jit(parallel=True)
def block_diag_loop_jit(arr, num=NUM):
rows, cols = arr.shape
result = np.zeros((num * rows, num * cols), dtype=arr.dtype)
for k in nb.prange(num):
for i in nb.prange(rows):
for j in nb.prange(cols):
result[i + (rows * k), j + (cols * k)] = arr[i, j]
return result针对NUM = 4的基准测试
针对NUM = 400的基准测试
绘图是使用以下附加代码从this template生成的：
def gen_input(n):
return np.random.randint(1, M, (n, n))
def equal_output(a, b):
return np.all(a == b)
funcs = block_diag_kron, block_diag_scipy, block_diag_view, block_diag_jit
input_sizes = tuple(int(2 ** (2 + (3 * i) / 4)) for i in range(13))
print('Input Sizes:\n', input_sizes, '\n')
runtimes, input_sizes, labels, results = benchmark(
funcs, gen_input=gen_input, equal_output=equal_output,
input_sizes=input_sizes)
plot_benchmarks(runtimes, input_sizes, labels, units='ms')(编辑以包括np.einsum()-based方法和另一个具有显式循环的Numba版本。)
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/33508322
复制相关文章
相似问题
腾讯云开发者
