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.

    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.

    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 ?

  • $\begingroup$ 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... $\endgroup$ Commented Apr 19, 2023 at 15:14
  • 1
    $\begingroup$ 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". $\endgroup$ Commented Apr 19, 2023 at 15:19
  • 1
    $\begingroup$ "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.) $\endgroup$ Commented Apr 19, 2023 at 16:21
  • 1
    $\begingroup$ and no, by bootstrapping, I mean resampling, computing an MLE (and no hessian), and repeating. Then, calculating empirical standard error on the bootstrapped sample. $\endgroup$ Commented Apr 19, 2023 at 16:22
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Apr 20, 2023 at 16:00

1 Answer 1


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.

Wikipedia cautions about this issue, but I think they understate it:

In statistical estimation problems (such as maximum likelihood or Bayesian inference), credible intervals or confidence intervals for the solution can be estimated from the inverse of the final Hessian matrix[citation needed]. However, these quantities are technically defined by the true Hessian matrix, and the BFGS approximation may not converge to the true Hessian matrix

  • $\begingroup$ Thank you for your answer. If I understand in short BFGS performs a rough and inconsistent approximation of the matrix that doesn't necessarily converges toward the real one. So one should use statmodel approx_hessian instead after a MLE to retrieve standard error I guess ? $\endgroup$ Commented Apr 19, 2023 at 15:31
  • 1
    $\begingroup$ Yep. As a small caveat, usually when we use BFGS we have a closed form solution for the gradient, but maybe not the hessian. It's usually a faster to estimate the hessian off first order differences of the gradient than to estimate it from second order differences of the target function as is done in approx_hessian, but I would only worry about this if speed is a significant issue and the parameter space is moderately sized (i.e. >10). $\endgroup$
    – Cliff AB
    Commented Apr 19, 2023 at 15:38
  • $\begingroup$ All right in my case computing the hessian itself is not a big deal (1 minute for 1500 observations in a 12 parameters space) but then the problem is : how do I deal with negative element of the hessian inverse diagonal ? There are not variance so I can't retrieve any standard error from them ? $\endgroup$ Commented Apr 19, 2023 at 15:44
  • 1
    $\begingroup$ If the hessian is not invertible, it means the likelihood is not well behaved and unfortunately MLE theory cannot be used to derive standard errors. Sometimes this issue will go away simply by having more data if that's available. In other cases, you can try to reparameterize the model so that the likelihood behaves better, i.e. reparameterize so that derivatives are more stable. $\endgroup$
    – Cliff AB
    Commented Apr 19, 2023 at 16:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.