我正在尝试实现kronecker product函数。以下是我有两个想法:
def kron(arr1, arr2):
"""columnwise outer product, avoiding relocate elements.
"""
r1, c1 = arr1.shape
r2, c2 = arr2.shape
nrows, ncols = r1 * r2, c1 * c2
res = np.empty((nrows, ncols))
for idx1 in range(c1):
for idx2 in range(c2):
new_c = idx1 * c2 + idx2
temp = np.zeros((r2, r1))
temp_kron = scipy.linalg.blas.dger(
alpha=1.0, x=arr2[:, idx2], y=arr1[:, idx1], incx=1, incy=1,
a=temp)
res[:, new_c] = np.ravel(temp_kron, order='F')
return res
def kron2(arr1, arr2):
"""First outer product, then rearrange items.
"""
r1, c1 = arr1.shape
r2, c2 = arr2.shape
nrows, ncols = r1 * r2, c1 * c2
tmp = np.outer(arr2, arr1)
res = np.empty((nrows, ncols))
for idx in range(arr1.size):
for offset in range(c2):
orig = tmp[offset::c2, idx]
dest_coffset = idx % c1 * c2 + offset
dest_roffset = (idx // c1) * r2
res[dest_roffset:dest_roffset+r2, dest_coffset] = orig
return res
基于此stackoverflow帖子,我创建了一个MeasureTime装饰器。一个自然的基准将是与numpy.kron进行比较。以下是我的测试功能:
@MeasureTime
def test_np_kron(arr1, arr2, number=1000):
for _ in range(number):
np.kron(arr1, arr2)
return
@MeasureTime
def test_kron(arr1, arr2, number=1000):
for _ in range(number):
kron(arr1, arr2)
@MeasureTime
def test_kron2(arr1, arr2, number=1000):
for _ in range(number):
kron2(arr2, arr1)
证明Numpy的kron函数性能要好得多:
arr1 = np.array([[1,-4,7], [-2, 3, 3]], dtype=np.float64)
arr2 = np.array([[8, -9, -6, 5], [1, -3, -4, 7], [2, 8, -8, -3], [1, 2, -5, -1]], dtype=np.float64)
我想知道为什么会这样吗?那是因为Numpy的重塑方法比手动复制内容(尽管仍然使用numpy)性能更高?我很困惑,因为否则我也在使用np.outer / blas.dger。我在这里认识到的唯一区别是我们如何安排最终结果。NumPy的重塑效果如何?
Here是NumPy 1.17 kron源的链接。
更新:首先忘了提到我试图在python中进行原型设计,然后使用带有cblas / lapack的C ++实现kron。有一些现有的“克朗”需要重构。然后,我遇到了Numpy的reshape并给人留下了深刻的印象。
提前感谢您的时间!
In [124]: A, B = np.array([[1,2],[3,4]]), np.array([[10,11],[12,13]])
kron
产生:
In [125]: np.kron(A,B)
Out[125]:
array([[10, 11, 20, 22],
[12, 13, 24, 26],
[30, 33, 40, 44],
[36, 39, 48, 52]])
outer
产生相同的数字,但排列不同:
In [126]: np.outer(A,B)
Out[126]:
array([[10, 11, 12, 13],
[20, 22, 24, 26],
[30, 33, 36, 39],
[40, 44, 48, 52]])
[kron
将其重塑为A
和B
的形状的组合:
In [127]: np.outer(A,B).reshape(2,2,2,2)
Out[127]:
array([[[[10, 11],
[12, 13]],
[[20, 22],
[24, 26]]],
[[[30, 33],
[36, 39]],
[[40, 44],
[48, 52]]]])
然后将concatenate
将4维重新组合为2:
In [128]: np.concatenate(np.concatenate(_127, 1),1)
Out[128]:
array([[10, 11, 20, 22],
[12, 13, 24, 26],
[30, 33, 40, 44],
[36, 39, 48, 52]])
一种替代方法是交换轴,然后重塑形状:
In [129]: _127.transpose(0,2,1,3).reshape(4,4)
Out[129]:
array([[10, 11, 20, 22],
[12, 13, 24, 26],
[30, 33, 40, 44],
[36, 39, 48, 52]])
第一次重塑和转置会产生一个视图,但是第二次重塑必须产生一个副本。串联会复制。但是所有这些操作都是在已编译的numpy
代码中完成的。
定义功能:
def foo1(A,B): temp = np.outer(A,B) temp = temp.reshape(A.shape + B.shape) return np.concatenate(np.concatenate(temp, 1), 1) def foo2(A,B): temp = np.outer(A,B) nz = temp.shape temp = temp.reshape(A.shape + B.shape) return temp.transpose(0,2,1,3).reshape(nz)
测试:
In [141]: np.allclose(np.kron(A,B), foo1(A,B))
Out[141]: True
In [142]: np.allclose(np.kron(A,B), foo2(A,B))
Out[142]: True
时间:
In [143]: timeit np.kron(A,B)
42.4 µs ± 294 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [145]: timeit foo1(A,B)
26.3 µs ± 38.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [146]: timeit foo2(A,B)
13.8 µs ± 19.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
我的代码可能需要一些概括,但是它证明了这种方法的有效性。
===与您的
kron
:
In [150]: kron(A,B) Out[150]: array([[10., 11., 20., 22.], [12., 13., 24., 26.], [30., 33., 40., 44.], [36., 39., 48., 52.]]) In [151]: timeit kron(A,B) 55.3 µs ± 1.59 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)