Python (numpy) question: I have two 2D numpy arrays, A[i,j] and B[k,l], but the indexes are unrelated to each other (A and B won't even have the same dimensions in general). I want to be able to add/multiply these two together to get a 4D matrix: C[i,j,k,l] = A[i,j] + B[k,l] or C[i,j,k,l] = A[i,j]*B[k,l]. Obviously I can do this by looping through all four indices and running the statement "C[i,j,k,l] = A[i,j] + B[k,l]" or "C[i,j,k,l] = A[i,j]*B[k,l]" at each iteration of the loop, but is there an elegant way of doing this using a numpy function rather than loops? Some sort of addition function that adds an N-dimensional array to an M-dimensional array to give an (N+M)-dimensional array of all possible pairwise additions of elements? My loops are taking a while to run, and I'm hoping that using a "proper" numpy array manipulation function rather than loops would speed things up (and clean up my code a bit).