Efficiently computing element wise product of transition matrices (m*m) * (n*n) to give (mn*mn) matrix

Consider input matrices X and Y of shapes (m,m) and (n,n) respectively. As an output we need to give a (mn,mn) shape matrix such that it multiplies corresponding entries in the two matrices.
These two matrices X and Y represent transition matrices. A following example can be taken to illustrate the required output. Here, X is a 33 matrix and Y is a 22 matrix.

Matrix X
    x1  x2  x3    
x1|  a   b   c
x2|  d   e   f
x3|  g   h   i

Matrix Y
    y1  y2
y1|  j   k
y2|  l   m

Matrix Z (Output)
      x1y1  x1y2  x2y1  x2y2  x3y1  x3y2
x1y1|  aj    ak    bj    bk    cj    ck
x1y2|  al    am    bl    bm    cl    cm
x2y1|  dj    dk    ej    ek    fj    fk


Following is a non-vectorized function I have written for this task:

def transition_multiply(X,Y):
    for i in range(num_rows_X):     
        for j in range(num_rows_Y):         
            for x in X[i]:
                 for y in Y[j]:                 
    return out

import numpy
print transition_multiply(numpy.array(X),numpy.array(Y))

I do get the required output but realize that the non-vectorized version would be really slow. What would be the best way to vectorize this computation, using Numpy.

To those who are interested why this computation is needed. It is needed in making transition matrix of a Factorial Hidden Markov Model from constituent transition matrices.

Best answer

This is Kronecker product, see here in numpy documentation.