It’s not so much that grad takes in scalar values (it can work with tensors of arbitrary size) - the limiting factor is that the autograd module can only compute jacobian-vector products. In case of taking grad of the gradients it’s actually computing a hessian-vector product of the original function. So, to recover the full hessian you have to do these N grad calls, where each call will multiply it by a one-hot vector (let’s say with a 1 in ith place), and will recover the ith row of the actual hessian.
I’m not familiar with the method for diagonal approximation, so I’d need to read up to answer your second question.
Hope this helps! Let me know if anything was unclear or you need any more help!