Solved by tom in post #2 
              
              
                You can use broadcasting. For example: 
y = x[:, :, None, :] * x[:, None, :, :]

There also is torch.einsum 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 (th…