Thanks for your thoughts @albanD ! So the implementation above is not the vanilla 2-backward calls implementation of HVPs, it’s supposed to implement the Pearlmutter trick, which is indeed mixing forward and backward AD if I understand correctly. We can replace the forward-mode with two reverse-mode calls (see for reference j-town’s  A new trick for calculating Jacobian vector products 60) so altogether three backward calls as above should be fine. Please correct me if I’m wrong.
Now what I don’t understand is why this approach is taking so much longer than both the naive (“Simple AD” above) approach and the finite different approximation. It’s one order of magnitude slower than finite difference approximation which requires two gradient evaluations and asymptotically speaking should have the same complexity as Pearlmutters method. I guess I’m lacking some knowledge on autograds internals here… 
EDIT: So rereading the above linked blog post’s conclusion actually brings some clarification:

For Autograd, the situation isn’t quite as simple. Autograd doesn’t do its differentiation ahead of time, and needs to trace a function’s execution in order to reverse-differentiate it. This makes our new trick, in its basic form, significantly less efficient than the forward mode implementation 27 that I’ve written.

So this is also the case for PyTorch and PyTorch’s autograd I guess? As a result the trick doesn’t bring any advantage because it computes more than “necessary” in some way? An actual forward-mode implementation is the only way to move forward?