I have two Numpy arrays
x
(m, i)
y
(m, j)
x
y
(m, i*j)
import numpy as np
np.random.seed(1)
x = np.random.randint(0, 2, (10, 3))
y = np.random.randint(0, 2, (10, 2))
x
array([[1, 1, 0],
[0, 1, 1],
[1, 1, 1],
[0, 0, 1],
[0, 1, 1],
[0, 0, 1],
[0, 0, 0],
[1, 0, 0],
[1, 0, 0],
[0, 1, 0]])
y
array([[0, 0],
[1, 1],
[1, 1],
[1, 0],
[0, 0],
[1, 1],
[1, 1],
[1, 1],
[0, 1],
[1, 0]])
array([[0, 0, 0, 0, 0, 0],
[0, 0, 1, 1, 1, 1],
[1, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 1, 1],
[0, 0, 0, 0, 0, 0],
[1, 1, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0]])
x
y
def _mult(x, y):
r = []
for xc in x.T:
for yc in y.T:
r.append(xc * yc)
return np.array(r).T
Use NumPy broadcasting

(y[:,None]*x[...,None]).reshape(x.shape[0],1)
Explanation
As inputs, we have 
y : 10 x 2
x : 10 x 3
With y[:,None]
, we are introducing a new axis between the existing two dims, thus creating a 3D
array version of it. This keeps the first axis as the first one in 3D
version and pushes out the second axis as the third one.
With x[...,None]
, we are introducing a new axis as the last one by pushing up the two existing dims as the first two dims to result in a 3D
array version.
To summarize, with the introduction of new axes, we have 
y : 10 x 1 x 2
x : 10 x 3 x 1
With y[:,None]*x[...,None]
, there would be broadcasting
for both y
and x
, resulting in an output array with a shape of (10,3,2)
. To get to the final output array of shape (10,6)
, we just need to merge the last two axes with that reshape.