Question: I have 2 numpy array, one is of shape `(2, 5, 10)` and the other one is of shape `(2, 5, 10, 10)`.
The multiplication I want is `[[x1[i][j] * x2[i][j] for j in range(5)] for i in range(2)]`
It works as intended but is pretty slow, and I wanted to multiply directly x1 \* x2 but numpy does not like that.
Is there a numpy method to multiply over a given axis?

I have tried the numpy.multiply for it says 'axis' is an invalid keyword to ufunc 'multiply'

```py
x1 = np.arange(100).reshape((2, 5, 10))
x2 = np.arange(1000).reshape((2, 5, 10, 10))
x = [[x1[i][j] * x2[i][j] for j in range(5)] for i in range(2)] # slow method that works
x = np.multiply(x1, x2, axis=2) # What I'm looking for but doesn't work.

```

Answer: 
One way of doing this, assuming that you are fine with `x` as type `np.array`, is by making better use of [NumPy broadcasting](https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html) (see @hpaulj's answer):

```py
import numpy as np

x1 = np.arange((2 * 3 * 4)).reshape((2, 3, 4))
x2 = np.arange((2 * 3 * 4 * 4)).reshape((2, 3, 4, 4))
x = [[x1[i][j] * x2[i][j] for j in range(3)] for i in range(2)]) # slow method that works

# : using NumPy broadcasting
y = x1[:, :, None, :] * x2

np.all(np.array(x) == y)
# True

```

Timewise is a ~10-fold acceleration:

```py
%timeit np.array([[x1[i][j] * x2[i][j] for j in range(3)] for i in range(2)])
# 13.6 µs ± 388 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit x1[:, :, None, :] * x2
# 1.4 µs ± 10.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

```
Using your input we have

`x1.shape
(2, 5, 10)
x2.shape
(2, 5, 10, 10)`

You can use np.tensordot to multiply along arbitrary axes - as long as they match in dimension.

`x3 = np.tensordot(x1, x2, (1,1))
x3.shape
(2, 10, 2, 10, 10)`

tensordot also allows you to multiply and add over a number of axes. In this case for example.

`x4 = np.tensordot(x1,x2, ((0,1), (0,1)))
x4.shape
(10, 10, 10)`

Here the argument ((0,1), (0,1)) refers to the axes of x1 and x2, respectively.