Hi,
that is a great observation and good research to pinpoint the exact line that is problematic!
I think this has two parts:

There seems to be indeed something going on with matmul that shouldn’t. The timing below suggests that matmul has 4.5ms while even the naive mul + sum has 0.6ms with the below:

     x = torch.rand(1000, 10, 3, 4)
     y = torch.rand(10, 4, 5)
     y2 = y.unsqueeze(0)
     x_ = x.unsqueeze(-1)
     y_ = y2.unsqueeze(2)
     %timeit z = torch.matmul(x, y)
     %timeit z2 = (x_*y_).sum(-2)


The strategy matmul employs is to reshape the array into something to be fed to at::bmm. Actually, in your example, it is not possible to do this, as you cannot combine the stride 0 of the first axis and the stride of the second axis into one. So a proper fix probably boils down to reworking the matmul rather fundamentally. (Probably a fun project!)

If you wanted to, you could file a bug. I do believe there is a bug in the code, not just the docs.
Best regards
Thomas