Python scipy.optimize.minimize: ‘trust-constr’ and Hessian output

Question

I'm using scipy.optimize.minimize, method = 'trust-constr'.

Does anyone know how to retrieve the Hessian at the minimum when using method = ‘trust-constr’? I would like to use the Hessian to compute standard errors of my estimates.

Answer 1

Lets assume you have formulated your problem (in this MCVE a least square problem):

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

You have some data and uncertainty on it:

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

When you minimize, it recovers the initial parameters:

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

What you need is to assess the Hessian of the function at the solution:

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

And invert the Hessian matrix to estimate the Covariance matrix:

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

From this you can draw uncertainty on each parameter:

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