Thank you again for your comment.

I’m not surprised that you only get 3 digits of precision (remember that float is 6 digits precision for a single operation)

What did you mean here? 3 digits of precision come from that, in the above example, I have
tensor([-0.7946])
tensor([-0.7947])
tensor([-0.7947], grad_fn=<SubBackward0>)

Is it correct?
So you might suspect that the difference between the automatic differentiation and the numerical approximation with torch.autograd.gradcheck() comes from a numerical precision issue. You might recommend me to use the Jacobian computed from torch.autograd.grad and use it in an optimization?