The “per-example gradient” computation needs to distinguish between “data” and “parameter” inputs, so it’s natural to handle them at the level of layers (ie, subclasses of nn.Module) rather than at the level of ops.
A practical limitation is that PyTorch ops often don’t support batch dimension for the weight input. For example, https://arxiv.org/abs/1906.02506 11 needed to compute loss over a batch of weights in PyTorch. The solution had to use 8 copies of model for a batch of size 8 (one per GPU).