Hi,
I am working on the same problem, I used the code of the link: https://gist.github.com/apaszke/226abdf867c4e9d6698bd198f3b45fb7 10, which you posted earlier at this post.
I am also using a straight forward neuronal network as a toy example:
def jacobian(y, x, create_graph=True):                                                               
    jac = []                                                                                          
    flat_y = y.reshape(-1)
    #print("flat_y: ", flat_y)  

    grad_y = torch.zeros_like(flat_y)   
    #print("grad_y: ", grad_y)    

    for i in range(len(flat_y)):                                                                      
        grad_y[i] = 1.                                                                                
        grad_x, = torch.autograd.grad(flat_y, x, grad_y, retain_graph=True, create_graph=create_graph)
        #grad_x, = torch.autograd.grad(flat_y, x,retain_graph=True, create_graph=create_graph)
        jac.append(grad_x.reshape(x.shape))                                                           
        grad_y[i] = 0.                                                                                
    return torch.stack(jac).reshape(y.shape + x.shape)                                                
   
                                                                                                      
def hessian(y, x):                                                                                    
    return jacobian(jacobian(y, x, create_graph=True), x) 

if __name__ == '__main__':
    x = torch.ones(4, requires_grad=True) 
    y = Variable(Tensor([1.]), requires_grad=True)

    model = nn.Linear(4, 1)
    model.weight.data.fill_(0.5)
    model.bias.data.fill_(0.5)

    loss = torch.sum((y - model(x))**2)
    pred = model(x)
 
    print("type", loss)
    print("type(loss)", type(loss))
    print("loss.shape", loss.shape)
    print(hessian(loss, x))
    print("\n") 

    print("pred", pred[0])
    print("type(pred)", type(pred[0]))
    print("pred.shape", pred[0].shape)
    print("\n")   
    print(hessian(pred[0], x))     



And get the following result and error:
loss:  tensor(2.2500, grad_fn=< SumBackward0 >)
type(loss):  <class ‘torch.Tensor’>
loss.shape:  torch.Size([])
tensor([[0.5000, 0.5000, 0.5000, 0.5000],
[0.5000, 0.5000, 0.5000, 0.5000],
[0.5000, 0.5000, 0.5000, 0.5000],
[0.5000, 0.5000, 0.5000, 0.5000]], grad_fn=< ViewBackward >)
pred:  tensor(2.5000, grad_fn=< SelectBackward >)
type(pred):  <class ‘torch.Tensor’>
pred.shape:  torch.Size([])
Traceback (most recent call last):
File “test_hfm.py”, line 58, in 
print(hessian(pred[0], x))
File “test_hfm.py”, line 27, in hessian
return jacobian(jacobian(y, x, create_graph=True), x)
File “test_hfm.py”, line 19, in jacobian
grad_x, = torch.autograd.grad(flat_y, x, grad_y, retain_graph=True, create_graph=create_graph)
File “/Users/felixasanger/anaconda3/envs/dil_bmw/lib/python3.7/site-packages/torch/autograd/init.py”, line 158, in grad
inputs, allow_unused)
RuntimeError: One of the differentiated Tensors appears to not have been used in the graph. Set allow_unused=True if this is the desired behavior.
Why does everything work fine for the loss, but not if I try it with the prediction? Is there any difference?
Any help is highly apperciated, I read a lot about this problem, but I just can´t solve it.
Best regards and thanks a lot,
Felix