If
a.shape = (r0,r1,..,rN) and b.shape = (s0,s1,...,sN), the Kronecker product has shape (r0*s0, r1*s1, ..., rN*SN).np.einsum('ij,kl->ikjl', A, B).reshape(A.shape[0] * B.shape[0], A.shape[1] * B.shape[1])numpy.kron — NumPy v2.1 Manual
Computes the Kronecker product, a composite array made of blocks of the
second array scaled by the first.
https://numpy.org/doc/stable/reference/generated/numpy.kron.html

Seonglae Cho