Let’s say f: R^N -> R, and f is batched so that f(xs) -> ys, where xs.shape == (n_batch, n_dim) and ys.shape = (n_batch,). I want to calculate the Laplacian of f, that is, \sum_i d^2/dx_i^2 f(x).
At the moment, I evaluate ys = f(xs), and then iterate over n_dim, and for each i a I calculate
diys = grad(ys, xis, ones(n_batch), create_graph=True)
didiys = grad(diys, xis, ones(n_batch), retain_graph=True)

and then sum all didiys. Is there a way to do this in on step rather than iterating over dimensions?
It is somewhat related to Issue #7786 16 where f_1, f_2  would be d/dx1 f(x), d/dx2 f(x).