Thank you for your answer!
Then, it appears that .backward() always “returns” a vector with the same dimension of the parameters (to be more specific, it saves in the .grad attribute of each parameter only one number. Either the derivative of a scalar loss with respect to this parameter, either the scalar product between a fixed vector vec (given as input in .backward()), and the vector made of all the derivatives with respect to this parameter of a vectorial loss).
In particular, there is no way to parallelize the computation of the Jacobian? I guess there’s not…
Furthermore, let’s say my vec I give as input in .backward() is just a vector of ones, since I’m interested in the sum of the gradients. Is it better to approach the problem this way, or to sum up the outputs of my vectorial loss before? (This should be equivalent for linearity, but maybe one method is computationally / numerically better…