This is a case for broadcasting: You add size-1 dimensions to gamma so it has ones for the last two axes: x * gamma[:, :, None, None] the : takes a full dimension from the original gamma, the None adds a singleton dimension which is then treated with broadcasting.
You could also pretend to be fancy by using einsum: torch.einsum('ijkl,ij->ijkl', x, gamma) will give the same result. If you remove letters from the last bit (following the ->) of the equation string, einsum will sum over these.