I have almost the same issue as you had. I have a model with input N=15 and output M=55, and I need to compute the Jacobian with respect to the input parameters (15).
Using torch.autograd.functional.jacobian is almost 8 times slower than the finite difference method.
I am curious whether autograd can be implemented in a way which is faster than finite difference.