I want to compute 3rd order partial derivatives of a scalar function w.r.t. an input matrix. So if the input matrix X is 5x5 (so 25 entries), I think I should end up with (5x5)x(5x5)x(5x5) entries (25^3), right?
Main question: Is this a correct and if so, the right/best/fastest/… way to do it?
def jacobian(func, inp):
    return torch.autograd.functional.jacobian(func, inp, create_graph=True)

def third_partial_derivatives(scalar_func, inp):
    return jacobian(
        lambda x: jacobian(
            lambda y: jacobian(scalar_func, y), 
            x
        ), 
        inp
    )

Side question: I noticed that the 2nd order derivatives already contain almost exclusively very small values, and the 3rd order derivatives contain almost only zeroes. Could this be caused by underflow, and is there a way I could fix this?
Thanks in advance!