This is because torch.autograd.functional.jacobian differentiates with respect to the whole Tensor, i.e. it includes the batch dim as a differentiable dimension of the input array (hence why you get the batch dim twice). You can then just take the diagonal over the pair of batch dimensions to get the per-sample grads you want! This can be done via torch.einsum
grads = torch.einsum("bibj->bij",grads)