Python scipy.optimize.minimize：‘trust-constr’与Hessian输出

Question

我正在使用scipy.optimize.minimize这个工具，方法是'trust-constr'。

有没有人知道在使用'trust-constr'这个方法时，怎么获取到最小值处的Hessian矩阵？我想用Hessian矩阵来计算我估计值的标准误差。

Answer 1

假设你已经把问题整理清楚了（在这个简单的例子中是一个最小二乘问题）：

import inspect
import numpy as np
import numdifftools as nd
from scipy import optimize

def model(x, a, b):
    return a * x + b

def loss_factory(func, x, y, s):
    ddof = len(inspect.signature(func).parameters) - 1
    def wrapped(p):
        return np.sum(np.power((y - func(x, *p)) / s, 2)) / (x.size - ddof)
    return wrapped

你有一些数据，并且对这些数据有一定的不确定性：

np.random.seed(123456)
x = np.linspace(-1, 1, 30)
s = 0.1 * np.ones_like(x)
y = model(x, 2, 3)
n = s * np.random.normal(size=x.size)
y += n

当你进行最小化操作时，它会恢复出最初的参数：

loss = loss_factory(model, x, y, s)
sol = optimize.minimize(loss, x0=[1., 1.], method="trust-constr")
# x: array([1.98896756, 2.96343096])   

你需要做的是在解决方案上评估这个函数的Hessian：

H = nd.Hessian(loss)(sol.x)
# array([[ 7.63546798e+01, -1.22821246e-14],
#        [-1.22821246e-14,  2.14285714e+02]])

然后反转Hessian矩阵，以估算协方差矩阵：

C = np.linalg.inv(H)
# array([[1.30967742e-02, 7.50662328e-19],
#        [7.50662328e-19, 4.66666667e-03]])

通过这个，你可以得出每个参数的不确定性：

sp = np.sqrt(np.diag(C))
# array([0.11444114, 0.06831301])