矩阵阵列的克罗内克积答案

【问题标题】：Kronecker product of matrix array矩阵阵列的克罗内克积
【发布时间】：2019-12-07 03:09:09
【问题描述】：

我有两个形状相同的矩阵数组 A 和 B： A.shape = B.shape = (M,N,P)

我想沿轴 0 计算克罗内克积，这样：

KP[ii,:,:] = A[ii,:,:]⊗B[ii,:,:]

有没有办法在 numpy 中做到这一点而不使用 for 循环？

谢谢！

例子：

A = np.array([ [[1,0],
                [0,1]],
               [[1,0],
                [0,1]]
            ])

B = np.array([ [[1,0],
                [0,-1]],
               [[0,1],
                [1,0]]
            ])

KP = np.array( [
                [[1,0,0,0],
                 [0,-1,0,0],
                 [0,0,1,0],
                 [0,0,0,-1]],
                [[0,1,0,0],
                 [1,0,0,0],
                 [0,0,0,1],
                 [0,0,1,0]]
               ] )

相当于：

KP= np.zeros( (A.shape[0],
               A.shape[1]**2,
               A.shape[2]**2) )

for ii in range(A.shape[0]):

    KP[ii,:,:] = np.kron(A[ii,:,:],B[ii,:,:])

【问题讨论】：

那么KP会有什么形状呢？
你能用(2, 2, 2) 数组之类的东西创建一个愚蠢的例子吗，我不确定你到底想要什么。
KP.shape = (M,N^2,P^2)
np.kron(A, B)[:A.shape[0]] - 这只是因为 A 由对称块组成
@SoslanTabuev 你的意思是相同的块？

标签： python-3.x numpy numpy-ndarray array-broadcasting

【解决方案1】：

你可以使用einsum，稍加练习就可以很直观，或者走经典的重塑和广播路线

A = np.array([ [[1,0],
                [0,1]],
               [[1,0],
                [0,1]]
            ])

B = np.array([ [[1,0],
                [0,-1]],
               [[0,1],
                [1,0]]
            ])

i,j,k = A.shape
i,l,m = B.shape
np.einsum("ijk,ilm->ijlkm",A,B).reshape(i,j*l,k*m)

# array([[[ 1,  0,  0,  0],
#         [ 0, -1,  0,  0],
#         [ 0,  0,  1,  0],
#         [ 0,  0,  0, -1]],
# 
#        [[ 0,  1,  0,  0],
#         [ 1,  0,  0,  0],
#         [ 0,  0,  0,  1],
#         [ 0,  0,  1,  0]]])

等价的非einsum表达式：

(A[:,:,None,:,None]*B[:,None,:,None,:]).reshape(i,j*l,k*m)

【讨论】：

谢谢，这绝对有效！不过，我对非 einsum 表达式在这种情况下的工作方式有点困惑。
@blackbody 我在 B 索引中添加了尾随冒号。它是多余的，但也许使它更具可读性？

【解决方案2】：

以防有人发现这也有用（但也许效率不高）。

def give_kr_prod(matrices):
    #matrices list of 2 (or more in principle) matrices
    while len(matrices) != 1:
        sm, smf=[],[]
        for ind in range(len(matrices)):
            sm.append(matrices[ind])
            if ind%2==1 and ind>0:
                smf.append(np.kron(*sm))
                sm=[]
        matrices = smf
    return matrices[0]

matrices = np.random.randn(8,2,2)

【讨论】：