Thank you for the answer.
If the autograd tools can only do Jacobian Vector products, then, in my opinion, it’s quite confusing that you are able to specify a matrix with shape [n,m] for the grad_outputs parameter when the output is a matrix. For example, in the above scenario if I do
grad = torch.autograd.grad(
                outputs=output,
                inputs=input,
                grad_outputs=torch.ones(output.size()),
                retain_graph=True, 
                create_graph=True
            )[0]

then the grad has now shape [BATCH_SIZE, DIMENSIONALITY].