具有单位矩阵和正则矩阵的高效 Kronecker 乘积 - NumPy/Python答案

【问题标题】：Efficient Kronecker product with identity matrix and regular matrix - NumPy/ Python具有单位矩阵和正则矩阵的高效 Kronecker 乘积 - NumPy/Python
【发布时间】：2017-11-11 16:47:44
【问题描述】：

我正在开发一个 python 项目并使用 numpy.我经常需要通过单位矩阵计算矩阵的 Kronecker 积。这些是我的代码中的一个很大的瓶颈，所以我想优化它们。我必须服用两种产品。第一个是：

np.kron(np.eye(N), A)

只需使用scipy.linalg.block_diag 就可以很容易地优化这个。该产品相当于：

la.block_diag(*[A]*N)

这大约快 10 倍。但是，我不确定如何优化第二种产品：

np.kron(A, np.eye(N))

我可以使用类似的技巧吗？

【问题讨论】：

标签： python numpy matrix scipy linear-algebra

【解决方案1】：

一种方法是初始化4D 的输出数组，然后从A 为其赋值。这样的赋值会广播值，这就是我们在 NumPy 中获得效率的地方。

因此，解决方案是这样的 -

# Get shape of A
m,n = A.shape

# Initialize output array as 4D
out = np.zeros((m,N,n,N))

# Get range array for indexing into the second and fourth axes 
r = np.arange(N)

# Index into the second and fourth axes and selecting all elements along
# the rest to assign values from A. The values are broadcasted.
out[:,r,:,r] = A

# Finally reshape back to 2D
out.shape = (m*N,n*N)

作为一个函数 -

def kron_A_N(A, N):  # Simulates np.kron(A, np.eye(N))
    m,n = A.shape
    out = np.zeros((m,N,n,N),dtype=A.dtype)
    r = np.arange(N)
    out[:,r,:,r] = A
    out.shape = (m*N,n*N)
    return out

要模拟np.kron(np.eye(N), A)，只需交换第一和第二轴的操作，第三和第四轴类似 -

def kron_N_A(A, N):  # Simulates np.kron(np.eye(N), A)
    m,n = A.shape
    out = np.zeros((N,m,N,n),dtype=A.dtype)
    r = np.arange(N)
    out[r,:,r,:] = A
    out.shape = (m*N,n*N)
    return out

时间安排 -

In [174]: N = 100
     ...: A = np.random.rand(100,100)
     ...: 

In [175]: np.allclose(np.kron(A, np.eye(N)), kron_A_N(A,N))
Out[175]: True

In [176]: %timeit np.kron(A, np.eye(N))
1 loops, best of 3: 458 ms per loop

In [177]: %timeit kron_A_N(A, N)
10 loops, best of 3: 58.4 ms per loop

In [178]: 458/58.4
Out[178]: 7.842465753424658

【讨论】：

啊，当然谢谢！作为参考，这在我的测试中提供了 6-7 倍的加速。
@user3930598 很高兴知道！添加了我的。似乎也在那里。