仅在 3 维数组的两个维度上的克罗内克积
Kronnecker product on just two dimensions of 3 dimensional arrays
假设有两个矩阵：
a = np.array([[1,2],
[3,4]])
b = np.array([[5,6],
[7,8]])
这将给出 Kronecker 乘积 ab = np.kron(a,b):
array([[ 1, 2, 2, 4],
[ 3, 4, 6, 8],
[ 3, 6, 4, 8],
[ 9, 12, 12, 16]])
现在假设在两个数组中有这些矩阵的三个副本，如下所示：
c = np.stack([a,a,a])
d = np.stack([b,b,b])
我想计算 c 和 d 的克罗内克积，这样输出是一个 3 索引数组，对应于 ab 的 3 个副本，即形状 [=17] =].然而，简单地执行 kron(c,d) 输出形状 (9,4,4)，其条目比需要的多，无法适当地重新整形。您能否帮助了解如何执行此操作？
假设您对中间列表没有问题（numpy 在分配数组之前需要大小），这是有效的：
import numpy as np
a = np.array([[1, 2],
[3, 4]])
b = np.array([[5, 6],
[7, 8]])
c = np.stack([a, a, a])
d = np.stack([b, b, b])
result = np.array(list(np.kron(x, y) for x, y in zip(c, d)))
print(result)
print(result.shape)
输出：
[[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]
[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]
[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]]
(3, 4, 4)
res=np.zeros((3,4,4))
res[:] = np.kron(a,b)
应该可以，向所有 3 架飞机广播 kron。
kron 是 a 和 b 外积的特殊重排。 A (2,2,2,2) 重新排列为 (4,4)。我在另一个 post:
中算出了细节
您的 (3,4,4) 可以从 (3,2,2,2,2) 获得，但它不是标准的，因此没有开箱即用的功能。您可以尝试修改我的答案。
In [246]: a = np.array([[1,2],
...: [3,4]])
...: b = np.array([[5,6],
...: [7,8]])
In [249]: np.kron(a,b)
Out[249]:
array([[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]])
正如我之前展示的，kron 可以通过对 outer 产品应用转置和形状来生成。我们可以对外部和转置使用 einsum：
In [253]: np.einsum('ij,kl->ikjl',a,b) # ikjl instead ijkl
Out[253]:
array([[[[ 5, 6],
[10, 12]],
[[ 7, 8],
[14, 16]]],
[[[15, 18],
[20, 24]],
[[21, 24],
[28, 32]]]])
In [254]: np.einsum('ij,kl->ikjl',a,b).reshape(4,4)
Out[254]:
array([[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]])
将其推广到数组 (3,2,2) 形状，我们可以添加一个额外的 'batch' 维度：
In [255]: c = np.stack([a,a,a])
...: d = np.stack([b,b,b])
In [256]: c
Out[256]:
array([[[1, 2],
[3, 4]],
[[1, 2],
[3, 4]],
[[1, 2],
[3, 4]]])
In [257]: c.shape
Out[257]: (3, 2, 2)
In [258]: np.einsum('aij,akl->aikjl',c,d).reshape(3,4,4)
Out[258]:
array([[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]],
[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]],
[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]]])
但如果我们知道 c 和 d 只是 a 和 b 的复制，则广播的解决方案更快
In [260]: res = np.zeros((3,4,4),int)
In [261]: res[:] = np.kron(a,b)
甚至更好（无副本的 kron 视图）：
np.broadcast_to(np.kron(a,b),(3,4,4))
一些时间：
In [280]: timeit np.einsum('aij,akl->aikjl',c,d).reshape(3,4,4)
10.2 µs ± 22.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [281]: timeit res=np.zeros((3,4,4),int);res[:] = np.kron(a,b)
47.5 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [282]: timeit np.broadcast_to(np.kron(a,b),(3,4,4))
57.6 µs ± 1.76 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [283]: timeit np.array(list(np.kron(x, y) for x, y in zip(c, d)))
143 µs ± 319 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
我有点惊讶 einsum 的速度有多快。 braodcast_to 也不是更快，尽管这主要是 kron 的错误（我之前的回答显示速度较慢）。
我的解决方案是这样的
import numpy as np
def kron(a, b):
ai = a[0]
bi = b[0]
b = np.asanyarray(b)
a = np.array(a, copy=False, subok=True, ndmin=b.ndim)
ndb, nda = bi.ndim, ai.ndim
b_size = a.shape[0]
as_ = ai.shape
bs = bi.shape
nd = ndb
if (ndb != nda):
if (ndb > nda):
as_ = (1,) * (ndb - nda) + as_
else:
bs = (1,) * (nda - ndb) + bs
nd = nda
res = np.einsum("ij,ik->ijk", a.reshape((a.shape[0], -1)), b.reshape((b.shape[0], -1)))
res = res.reshape((b_size, )+(as_ + bs))
axis = nd - 1
result = []
for i in range(b_size):
r = res[i]
for _ in range(nd):
r = np.concatenate(r, axis=axis)
result.append(r)
return np.array(result)
x = np.array([[1,2],
[3,4]])
y = np.array([[5,6],
[7,8]])
c = np.stack([x,x,x])
d = np.stack([y,y,y])
k = kron(c, d)
result = np.array(list(np.kron(x, y) for x, y in zip(c, d)))
print(np.allclose(k, result))
正如@hpaulj 在他们的回答中提到的那样，没有内置函数可以执行此操作，您必须对外部产品执行一些重要的重塑操作才能使其正常工作。
解决方案已根据 kron in numpy 的实现进行了修改。
我基本上只是在 numpy 实现用于执行串联操作的循环之上添加了另一个外部循环。这是为了分别处理 c 和 d 元素之间的每个外积。其他都差不多。
请注意，这不是 numpy 的 kron 实现的概括，因为这 不适用于 您只想直接执行 kron 的情况手术。至少不是没有一些额外的修改。为了清楚起见，我还忽略了 numpy 正在处理的一些边缘情况。但是，这应该适用于您的用例。
假设有两个矩阵：
a = np.array([[1,2],
[3,4]])
b = np.array([[5,6],
[7,8]])
这将给出 Kronecker 乘积 ab = np.kron(a,b):
array([[ 1, 2, 2, 4],
[ 3, 4, 6, 8],
[ 3, 6, 4, 8],
[ 9, 12, 12, 16]])
现在假设在两个数组中有这些矩阵的三个副本，如下所示：
c = np.stack([a,a,a])
d = np.stack([b,b,b])
我想计算 c 和 d 的克罗内克积，这样输出是一个 3 索引数组，对应于 ab 的 3 个副本，即形状 [=17] =].然而，简单地执行 kron(c,d) 输出形状 (9,4,4)，其条目比需要的多，无法适当地重新整形。您能否帮助了解如何执行此操作？
假设您对中间列表没有问题（numpy 在分配数组之前需要大小），这是有效的：
import numpy as np
a = np.array([[1, 2],
[3, 4]])
b = np.array([[5, 6],
[7, 8]])
c = np.stack([a, a, a])
d = np.stack([b, b, b])
result = np.array(list(np.kron(x, y) for x, y in zip(c, d)))
print(result)
print(result.shape)
输出：
[[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]
[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]
[[ 5 6 10 12]
[ 7 8 14 16]
[15 18 20 24]
[21 24 28 32]]]
(3, 4, 4)
res=np.zeros((3,4,4))
res[:] = np.kron(a,b)
应该可以，向所有 3 架飞机广播 kron。
kron 是 a 和 b 外积的特殊重排。 A (2,2,2,2) 重新排列为 (4,4)。我在另一个 post:
您的 (3,4,4) 可以从 (3,2,2,2,2) 获得，但它不是标准的，因此没有开箱即用的功能。您可以尝试修改我的答案。
In [246]: a = np.array([[1,2],
...: [3,4]])
...: b = np.array([[5,6],
...: [7,8]])
In [249]: np.kron(a,b)
Out[249]:
array([[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]])
正如我之前展示的，kron 可以通过对 outer 产品应用转置和形状来生成。我们可以对外部和转置使用 einsum：
In [253]: np.einsum('ij,kl->ikjl',a,b) # ikjl instead ijkl
Out[253]:
array([[[[ 5, 6],
[10, 12]],
[[ 7, 8],
[14, 16]]],
[[[15, 18],
[20, 24]],
[[21, 24],
[28, 32]]]])
In [254]: np.einsum('ij,kl->ikjl',a,b).reshape(4,4)
Out[254]:
array([[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]])
将其推广到数组 (3,2,2) 形状，我们可以添加一个额外的 'batch' 维度：
In [255]: c = np.stack([a,a,a])
...: d = np.stack([b,b,b])
In [256]: c
Out[256]:
array([[[1, 2],
[3, 4]],
[[1, 2],
[3, 4]],
[[1, 2],
[3, 4]]])
In [257]: c.shape
Out[257]: (3, 2, 2)
In [258]: np.einsum('aij,akl->aikjl',c,d).reshape(3,4,4)
Out[258]:
array([[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]],
[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]],
[[ 5, 6, 10, 12],
[ 7, 8, 14, 16],
[15, 18, 20, 24],
[21, 24, 28, 32]]])
但如果我们知道 c 和 d 只是 a 和 b 的复制，则广播的解决方案更快
In [260]: res = np.zeros((3,4,4),int)
In [261]: res[:] = np.kron(a,b)
甚至更好（无副本的 kron 视图）：
np.broadcast_to(np.kron(a,b),(3,4,4))
一些时间：
In [280]: timeit np.einsum('aij,akl->aikjl',c,d).reshape(3,4,4)
10.2 µs ± 22.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
In [281]: timeit res=np.zeros((3,4,4),int);res[:] = np.kron(a,b)
47.5 µs ± 1.66 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [282]: timeit np.broadcast_to(np.kron(a,b),(3,4,4))
57.6 µs ± 1.76 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [283]: timeit np.array(list(np.kron(x, y) for x, y in zip(c, d)))
143 µs ± 319 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
我有点惊讶 einsum 的速度有多快。 braodcast_to 也不是更快，尽管这主要是 kron 的错误（我之前的回答显示速度较慢）。
我的解决方案是这样的
import numpy as np
def kron(a, b):
ai = a[0]
bi = b[0]
b = np.asanyarray(b)
a = np.array(a, copy=False, subok=True, ndmin=b.ndim)
ndb, nda = bi.ndim, ai.ndim
b_size = a.shape[0]
as_ = ai.shape
bs = bi.shape
nd = ndb
if (ndb != nda):
if (ndb > nda):
as_ = (1,) * (ndb - nda) + as_
else:
bs = (1,) * (nda - ndb) + bs
nd = nda
res = np.einsum("ij,ik->ijk", a.reshape((a.shape[0], -1)), b.reshape((b.shape[0], -1)))
res = res.reshape((b_size, )+(as_ + bs))
axis = nd - 1
result = []
for i in range(b_size):
r = res[i]
for _ in range(nd):
r = np.concatenate(r, axis=axis)
result.append(r)
return np.array(result)
x = np.array([[1,2],
[3,4]])
y = np.array([[5,6],
[7,8]])
c = np.stack([x,x,x])
d = np.stack([y,y,y])
k = kron(c, d)
result = np.array(list(np.kron(x, y) for x, y in zip(c, d)))
print(np.allclose(k, result))
正如@hpaulj 在他们的回答中提到的那样，没有内置函数可以执行此操作，您必须对外部产品执行一些重要的重塑操作才能使其正常工作。
解决方案已根据 kron in numpy 的实现进行了修改。
我基本上只是在 numpy 实现用于执行串联操作的循环之上添加了另一个外部循环。这是为了分别处理 c 和 d 元素之间的每个外积。其他都差不多。
请注意，这不是 numpy 的 kron 实现的概括，因为这 不适用于 您只想直接执行 kron 的情况手术。至少不是没有一些额外的修改。为了清楚起见，我还忽略了 numpy 正在处理的一些边缘情况。但是，这应该适用于您的用例。