Reverse mode autograd (what we have in pytorch) is capable of computing vector-Jacobian products. That is, given a function f, an input x, and an arbitrary vector v, autograd can tell you v J where J is the Jacobian of f with x.
With multiple vector Jacobian products (If we query autograd for u1 J, u2 J, u3 J, etc where u is a unit vector) we can reconstruct the full Jacobian, which is what the code I wrote above does.