I was facing the same problem as OP and managed to get an elegant workaround using simple sum along the batch dimension of the original function. This means the result will be differentiated exactly as OP requested (just write it down). This approach is theoretically batch_size times faster than computing full jacobian and then using diagonal. I can see similar speedup in my calculations.
def batch_jacobian(f, x):
    f_sum = lambda x: torch.sum(f(x), axis=0)
    return jacobian(f_sum, x)

Note that it is necessary to permute the result to get the batch_dimension in front. In OPs case permute(1,0,2) will do the job.
The same trick can be applied for hessian and will bring even more significant speedup.