你可以完全按照你在最后一行说的做：对第二列的每一行进行单独的克罗内克积，然后将结果连接起来。
假设两个矩阵分别为x（4 x 3）和y（2 x 3），首先要做的是将x拆分为两部分，因为乘积的每一部分都只有一半矩阵参与。

x = x.reshape(2, 2, 3)

然后，您可以分别计算这两个产品：

z0 = np.kron(x[0], y[0])
z1 = np.kron(x[1], y[1])

最后，沿着第一个轴连接两个结果：

z = np.concatenate([z0, z1], axis=0)

或者，如果你像我一样，喜欢大丑俏皮话，你可以这样做：

z = np.concatenate([np.kron(xr, yr) for xr, yr in zip(x.reshape(2, 2, 3), y)], axis=0)

在你在评论中提到的一般情况下，它将变成：

z = np.concatenate([np.kron(xr, yr) for xr, yr in zip(x.reshape(int(n / 2), 2, 3), y)], axis=0)

这给显式循环带来了相同的结果，我相信它可以被numba.jit编译：

def solve_explicit(x, y):
    # sanity checks
    assert x.shape[0] == 2*y.shape[0]
    assert x.shape[1] == y.shape[1]

    n = x.shape[0]
    z = np.zeros((n, 9))
    for i in range(n):
        for j in range(3):
            for k in range(3):
                z[i, k + 3 * j] = x[i, j] * y[int(i / 2), k]
    return z

使用x.shape（n，3）和y.shape（n//2，3）的广播：

out = (x.reshape(-1, 2, 3, 1) * y.reshape(-1, 1, 1, 3)).reshape(-1, 9)

我个人在这种情况下会使用np.einsum，因为我认为它比广播更容易理解。

import numpy as np

(a, b, c, d) = np.random.rand(4)
x = np.array([[1, 0, a], [0, 1, b], [1, 0, c], [0, 1, d]])
y = np.random.rand(2, 3)
z = np.einsum("ij,ik->ijk", x.reshape(-1, 6), y).reshape(-1, 9)

# timeit magic commands.
# %timeit -n 50000 np.einsum("ij,ik->ijk", x.reshape(-1, 6), y).reshape(-1, 9)
# %timeit -n 50000 (x.reshape(-1, 2, 3, 1) * y.reshape(-1, 1, 1, 3)).reshape(-1, 9)

NumPy中关于爱因斯坦求和的一些很好的参考文献：[ 2，3，4 ]的第一个字符串。