Well, just the output matrix is going to be very expensive as it will take 4GB. So working with it will require a machine with a lot of memory.
But to go back to you original question, there is nothing much better than repetitive call to autograd.grad I’m afraid. As automatic differentiation only gives a row of the Jacobian at a time.