I would like to compute the dot product (matmul) of the inner dimension of two 3D arrays. In the following example, I have an array of 10 2x3 matrixes (X) and an array of 8 1x3 matrixes. The result Z should be a 10 element array of an 8 x 2 matrix (you might also think of this as an 10 x 8 array of 2-d vectors.) X = np.arange(0, 10 * 2 * 3).reshape(10, 2, 3) Y = np.arange(0, 8 * 1 * 3).reshape(8, 1, 3) # write to Z, which has the dot-product of internal dims of X and Y Z = np.empty((10, 8, 2)) for i, x_i in enumerate(X): for j, y_j in enumerate(Y): z = x_i @ y_j.T # subscriptz to flatten it. Z[i, j, :] = z[:, 0] is correct, but I would like a vectorized solution. Continue reading...