# scipy minimize gives a hess_inv that is completely different from inv(statmodel.approx_hess)

I'm fitting a model with MLE using scipy.minimize (method BFGS). I want to have the hessian to compute its inverse and retrieve the standard error of each parameter.

The result of minimize contains one approximation of the inverse hessian matrix in res.hess_inv. To cross-check with other methods, I computed the hessian inverse using statmodel approx_hess and then applying np.linalg.inv. It gives a completely different matrix (with some diagonal elements that have up to 2 order of magnitude difference)

Finally, I decided to compute the hessian myself since it looks like I can't cross-check the result of the functions above :

def approx_hessian(x, f, eps=1e-5, only_diag=False) -> np.ndarray:
"""
Computes an approximation of the Hessian matrix of a function using the finite difference method.

Parameters:
x (array-like): The point at which to compute the Hessian matrix.
f (callable): The function to compute the Hessian of.
The function should take an array-like object as input and return a scalar.
epsilon (float, optional): The step size for the finite difference method.

Returns:
hessian (ndarray): An approximation of the Hessian matrix of func at x.
"""
n = len(x)
hessian = np.zeros((n, n))

# Compute the off-diagonal elements
if not only_diag:
for i in range(n):
for j in range(n):
if i == j or i < j:
continue # the hessian is symmetrical, so we can compute only half of it
p1 = x.copy()
p1[i] += eps
p1[j] += eps

p2 = x.copy()
p2[i] += eps
p2[j] -= eps

p3 = x.copy()
p3[i] -= eps
p3[j] += eps

p4 = x.copy()
p4[i] -= eps
p4[j] -= eps

hessian[i][j] = (f(p1) - f(p2) - f(p3) + f(p4)) / (4 * eps ** 2)
hessian = hessian + hessian.transpose()  # fill the element under the diagonal with the same numbers

# compute diagonal elements
for i in range(n):
p_forward = x.copy()
p_forward[i] += eps

p_backward = x.copy()
p_backward[i] -= eps

hessian[i][i] = (f(p_forward) - 2 * f(x) + f(p_backward)) / (eps ** 2)

return hessian


I test my implementation by computing manually double derivatives on a simple function having 3 parameters :

def test_approx_hessian(self):
# example taken from https://www.allmath.com/hessian-matrix-calculator.php
def fn(params: np.ndarray) -> float:
x, y, z = params
return 2 * y * x ** 4 + 2 * y ** 3 + 2 * x * z ** 3

def expected_hessian(p: np.ndarray) -> np.ndarray:
x, y, z = p
return np.array([
[24 * x ** 2 * y, 8 * x ** 3, 6 * z ** 2],
[8 * x ** 3, 12 * y, 0],
[6 * z ** 2, 0, 12 * x * z],
])

point = np.array([1.0, 2.0, 3.0])
hessian_expected = expected_hessian(point)
hessian_approx = approx_hessian(point, fn)
self.assertTrue(np.allclose(hessian_expected, hessian_approx, atol=1e-4))


This test passes, so I believe my function works.

However when I compute the hessian inverse using my implementation, it gives a matrix that is not only still completely different from the 2 others, but also with diagonal elements that are negative ... (while they should be variances, so the diagonal elements are supposed to be positive as far as I know)

Are maths doing me a prank ? What am I doing wrong ? What function should I "trust" to have finally a standard error of my parameters ?

• couple things going on here: 1) you will only have positive variances if the Hessian is positive definite. If you have truly found on optimum, the hessian will be positive definite there. In high-dimensional problems, it is very easy to find approximate optima which don't actually have a positive definite hessian. in your final code block, it seems like you're evaluating the hessian at an arbitrary point. There is no reason to expect this to be positive definite. 2) BFGS only approximates the Hessian after many many runs, and only if the Hessian is positive definite at the stationary point... Commented Apr 19, 2023 at 15:14
• and finally, regarding " What function should I "trust" to have finally a standard error of my parameters ?" I would not trust any code not validated by simulation study (with known parameters where you can empirically evaluate the purported standard errors), and even then, we would still "trust but verify". Commented Apr 19, 2023 at 15:19
• "Likelihood is always concave at optimum right" yes, theoretically. But in practice, when optimizing nonconvex likelihoods in higher dimensional spaces using numerical methods, the stationary points we end up at tend not to be strict optima, but rather saddle points wrt at least some of the input directions. And yes that's right regarding the positive definite hessian; google "saddle point" to see visually what's going on. (This is all assuming that the eigenval is indeed negative, which we're not 100% sure about.) Commented Apr 19, 2023 at 16:21
• and no, by bootstrapping, I mean resampling, computing an MLE (and no hessian), and repeating. Then, calculating empirical standard error on the bootstrapped sample. Commented Apr 19, 2023 at 16:22
• I see. That makes things much more interesting :) I think solving this problem is nontrivial, and you may wish to seek an academic statistical collaborator (if you're in academia) or a statistical consultant (if you're in industry/gov't) rather than an answer on a site like CV. Commented Apr 20, 2023 at 16:00

You should not use the output from BFGS as the true inverse of the hessian.

The reason for this is that BFGS does not approximate the hessian based on the current rate of change of the derivatives of the target function. Rather, it builds up a useful approximation by tracking the change of the first derivative at each state of the optimization algorithm. I believe it also uses some regularization to keep the step sizes in check. Note that without regularization, you can't have an invertable Hessian based on solely first order differences of the gradient if you've got less steps than parameters of your model.

This means, at best, you have something like the inverse of the weighted average of the changes of the first derivative through the different iterations of your algorithm, which may be very different than the inverse of the Hessian at the MLE. At worst, you have a matrix that's heavily regularized as to make sure that your next step isn't too large, where the regularization has nothing to do with the uncertainty of the parameters of your model.