Thanks for the details.
I don’t have a full answer but the remarks I would have are:

These indexing ops do not get as much love in PyTorch as core NN code (like linear/conv) and so we might be able to optimize them. If you find one in particular that is especially lacking, you can open an issue and we should be able to improve that.
For small size, we know that PyTorch has much more “framework overhead” than numpy. This is mostly due to the fact that we have extra features (like autograd, different device backend, etc) that make going to the low level code a bit longer (order of microseconds). This is not visible for large ops, but very small ops in a loop will suffer heavily from that.