numpy中Kronecker产品的双循环矢量化和重塑答案

【问题标题】：Double loop vectorization and reshape for Kronecker product in numpynumpy中Kronecker产品的双循环矢量化和重塑
【发布时间】：2020-04-08 06:08:45
【问题描述】：

我有两个矩阵

import numpy as np
n = 10
a = 2*np.ones((n,n,3))
b = 3*np.ones((n,n,3))

我想以提醒克罗内克乘积的方式将它们相乘，然后总结

s = 0
for i in range(n):
    for j in range(n):
        s +=  a*b[i,j]

在 numpy 中是否存在对其进行矢量化的方法？

【问题讨论】：

您的代码正确吗？ a*b[i,j] 应该返回一个与a 形状相同的数组，而不是一个标量（正如您的s = 0 所建议的那样）。
是的，无论如何它都能正常工作，s 应该与a 的形状相同。但是也可以写s = np.zeros((n,n,3))

标签： python numpy matrix-multiplication tensor

【解决方案1】：

也许这可以用np.einsum()写得更优雅：

import numpy as np


n = 10
a = 2 * np.ones((n, n, 3))
b = 3 * np.ones((n, n, 3))

s = 0
for i in range(n):
    for j in range(n):
        s +=  a * b[i, j]

print(s.shape)
# (10, 10, 3)


ss = a * np.einsum('ijk->k', b)
print(ss.shape)
# (10, 10, 3)
print(np.all(s == ss))
# True

甚至只有np.sum():

sss = a * np.sum(b, axis=(0, 1))
print(sss.shape)
# (10, 10, 3)
print(np.all(s == sss))
# True

但np.einsum() 似乎更快：

n = 100
a = 2 * np.ones((n, n, 3))
b = 3 * np.ones((n, n, 3))

%timeit f_with_loops(a, b)
# 1 loop, best of 3: 787 ms per loop
%timeit a * np.einsum('ijk->k', b)
# 10000 loops, best of 3: 121 µs per loop
%timeit a * np.sum(b, axis=(0, 1))
# 1000 loops, best of 3: 254 µs per loop

【讨论】：

【解决方案2】：

你的代码可以改写为：

因此，这应该有效：

s = a * np.sum(np.sum(b,axis=1),axis=0)

【讨论】：

注意这里可以简单写成a * np.sum(b, axis=(0, 1))
谢谢。我什至不知道这个功能。