This should get you started:
import numpy as np
import itertools as it
def row_product(*arrays):
lengths = np.array([x.shape[0] for x in arrays])
positions = np.cumsum(lengths)
ranges = np.arange(positions[-1])
ranges = np.split(ranges,positions[:-1])
total = np.concatenate((arrays),axis=0)
inds = np.fromiter(it.chain.from_iterable(it.product(*ranges)), np.int)
inds = inds.reshape(-1, len(arrays))
return np.take(total, inds, axis=0)
The last dimension(s) must be the same.
Showing the results:
a=np.array([[2,0],[1,1],[2,0]])
b=np.array([[1,0],[0,1]])
c=np.array([[0,0]])
print row_product(a,b,c)
[[[2 0]
[1 0]
[0 0]]
[[2 0]
[0 1]
[0 0]]
[[1 1]
[1 0]
[0 0]]
[[1 1]
[0 1]
[0 0]]
[[2 0]
[1 0]
[0 0]]
[[2 0]
[0 1]
[0 0]]]
This is a 3D array where the unique combinations are in the last two axes. Seems to be reasonably fast, 1M unique combinations takes about 1/6 of a second.