Generalized Kronecker product with different type of product in numpy or scipy

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))
     )