Thank @AlbanOdot and @glaringlee for your answers. I will try using the autograd.grad 2 C++ function.
Another option I figured out is by following the C++ documentation, where is described how Autograd 1 works. However, the example reported there is not working since the variable c is not a scalar:
#include <torch/csrc/autograd/variable.h>
#include <torch/csrc/autograd/function.h>

torch::Tensor a = torch::ones({2, 2}, torch::requires_grad());
torch::Tensor b = torch::randn({2, 2});
auto c = a + b;
c.backward(); // a.grad() will now hold the gradient of c w.r.t. a.

In order to make this work, one has to call the backward function on each component,e.g.:
c[0][0].backward()

and then reset the gradients to obtain with a subsequent backward call the gradients w.r.t. a different output component:
a.grad().zero_();

I think the documentation should be updated to deal with the fact that backward on a non-scalar variable is not allowed anymore in LibTorch, since version 1.3 I think, as it was already in the Python version.
Best regards,
Luigi