Calculate Hessian matrix of a layer w.r.t loss for a batch of samples


I want to calculate the Hessian matrix for a batch of k inputs for a layer of a model. I know how to calculate it elementwise via nested gradient tapes, but so far and after reading the documentation I could not find how to vectorize it for a batch of input.

For example, if the batch is of size 13 and the layer is a (10,20) Dense layer, the output would be of shape (13,10,20,10,20).