Hi, I think to this day the only way is to use grad 337 function, but you will need to call it j times (once for each output). Unfortunately this requires many backward propagations and scales terribly with the output space of the function.
I have this exact same need, is there a different way to get the Jacobian of a function?