Slicing a big NumPy array
A one liner for your problem:

def slice_and_multiply(arr):
    return np.product(arr, axis=1, keepdims=True) / M


This multiplies all of the "columns" first, then divides by each "column" to get the result of multiplying all but that column, which may be problematic depending on your data due to potential overflows or precision losses. And it breaks completely if there are any zeros in the input array.

If you can get rid of the division (this is a typical interview question, by the way), then all of those problems dissappear. And you can do it with relative ease, by cleverly storing accumulated products from both ends of the array. then multiplying them together:

def slice_and_multiply_bis(arr):
    ret = np.ones_like(arr)
    np.cumprod(arr[:, :-1], axis=1, out=ret[:, 1:])
    ret[:, :-1] *= np.cumprod(arr[:, :0:-1], axis=1)[:, ::-1]
    return ret


And of course:

&gt;&gt;&gt; M = np.random.rand(16, 8, 1000, 255)
&gt;&gt;&gt; np.allclose(slice_and_multiply(M), slice_and_multiply_bis(M))
True


That last code is a little too clever for its own good, so it should it should be heavily commented.

By my timings, for your target size, slice_and_multiply is 5x faster than your implementation, and slice_and_multiply_bis about 2-3x faster.