from itertools import product
r, c = A.shape
arr_i = np.array(list(product(list(range(r)), repeat=n)))
arr_j = np.array(list(product(list(range(c)), repeat=n)))
out = A.ravel()[(arr_i * c)[:, (None), :] + arr_j].prod(2)
Numpy/Python: Efficient matrix as multiplication of cartesian product of input matrix