使用python,pytorch求海森Hessian矩陣

　　考慮一個函數$y=f(\textbf{x}) (R^n\rightarrow R)$，y的Hessian矩陣定義如下：

　　考慮一個函數:$$f(x)=b^Tx+\frac{1}{2}x^{T}Ax\\其中b^T=[1,3,5], A在代碼中可讀，可以自定義$$ 

　　求該函數在x = [0,0,0]處海森矩陣值的python代碼如下：

　　本代碼需要用到torch.autograd包中的核心函數torch.autograd.grad。相鄰隨筆中有詳細參考解析。大致原理是人工求導並保留了計算圖，所以求二階導很方便。如果想求解析形式的海森矩陣或許matlab中有相關的包能實現。

import torch


# 定義函數
x = torch.tensor([0., 0, 0], requires_grad=True)
b = torch.tensor([1., 3, 5])
A = torch.tensor([[-5, -3, -0.5], [-3, -2, 0], [-0.5, 0, -0.5]])
y = b@x + 0.5*x@A@x

# 計算一階導數,因為我們需要繼續計算二階導數,所以創建並保留計算圖
grad = torch.autograd.grad(y, x, retain_graph=True, create_graph=True)
# 定義Print數組,為輸出和進一步利用Hessian矩陣作准備
Print = torch.tensor([])
for anygrad in grad[0]:  # torch.autograd.grad返回的是元組
    Print = torch.cat((Print, torch.autograd.grad(anygrad, x, retain_graph=True)[0]))
print(Print.view(x.size()[0], -1))

　　輸出結果如下：