You can use broadcasting. For example:
y = x[:, :, None, :] * x[:, None, :, :]

There also is torch.einsum 57 for fancy stuff (torch.einsum("iaj,ibj->iabj", [x, x]), works in PyTorch self-compiled master), but in your case the above seems to be what looks most straightforward and transparent to me (the trailing “:” indices are for my taste, not technically necessary, alternatives include inserting the dimensions with unsqueeze instead of using indexing).
Best regards
Thomas