Consider two boolean arrays import numpy as np A = np.asarray([[True, False], [False, False]]) B = np.asarray([[False, True], [True, True]])  I want to take the kronecker product of A and B under the xor operation. The result should be: C = np.asarray([[True, False, False, True], [False, False, True, True], [False, True, False, True], [True, True, True, True]])  More generally, is there a simple way to implement the Kronecker product with some multiplication operator distinct from the operator *, in this instance the xor operator ^?
You could use broadcasting and reshaping: m, n = A.shape p, q = B.shape C = (A[:, None, :, None] ^ B[None, :, None, :]).reshape(m*p, n*q)  Simplified: C = (A[:, None, :, None] ^ B[None, :, None, :] ).reshape(A.shape[0]*B.shape[0], -1)  Also equivalent to: C = (np.logical_xor.outer(A, B) .swapaxes(1, 2) .reshape(A.shape[0]*B.shape[0], -1) )  Or with explicit alignment using repeat/tile without reshaping: p, q = B.shape C = np.repeat(np.repeat(A, p, axis=0), q, axis=1) ^ np.tile(B, A.shape)  Output: array([[ True, False, False, True], [False, False, True, True], [False, True, False, True], [ True, True, True, True]])  ND-generalization for N dimensional inputs, one could follow the same logic by expanding the dimensions in an interleaved fashion with expand_dims, before reshaping to the element-wise product of the dimensions: C = ( np.expand_dims(A, tuple(range(1, A.ndim*2, 2))) ^ np.expand_dims(B, tuple(range(0, A.ndim*2, 2))) ).reshape(np.multiply(A.shape, B.shape))  Interestingly, this is how kron is actually implemented in numpy (with some extra checks in place). Variant with outer: C = (np.logical_xor.outer(A, B) .transpose(np.arange(A.ndim+B.ndim) .reshape(-1, 2, order='F') .ravel()) .reshape(np.multiply(A.shape, B.shape)) ) 