Great solution @parthe - thanks! After testing it for correctness & speed I made these changes:
def get_sum_of_gradients(x):
    h = f(x)
    return torch.autograd.grad(h, x, create_graph=True)[0].sum(0)

Hessian_per_sample = jacobian(get_sum_of_gradients, x, vectorize=True).swapaxes(0, 1)

The alteration to get_sum_of_gradients was purely to avoid memory leak warnings.
The addition of vectorize=True led to significant speed-ups (~10x) for the cases I tested.
As an example test-case, I computed the Hessian of f(x), where f is a scalar-valued 3-layer feed forward neural network. I used a batch-size of 500 and the input / hidden dims of 200. Computation was on GPU. I got the following results:

grad + Jacobian (vectorize=True) : 0.02 seconds
grad + Jacobian (vectorize=False) : 0.22 seconds
For-loop over batch + Hessian (vectorize=True) : 7.24 seconds
For-loop over batch + Hessian (vectorize=False): 74.53 seconds