calculating the jacobian for a least-squares problem

I would like to calculate the Jacobian for a least squares problem,
followed by a Hessian estimation, then the covariance matrix from that
Hessian.

With my current approach I sometimes experience issues with the covariance
matrix in that it's sometimes not positive semi-definite. I am using the
covariance matrix to seed a MCMC sampling process by supplying it to
`np.random.multivariate_normal` to get initial positions for the MC chain.
I am using the following code:

```
from scipy.optimize._numdiff import approx_derivative
jac = approx_derivative(residuals_func, x0)
hess = np.matmul(jac.T, jac)
covar = np.linalg.inv(hess)
```

Note that x0 may not be at a minimum.

- would this be the usual way of estimating the Hessian, is there anything
incorrect with the approach?
- what is the recommended way (i.e. numerically stable) of inverting the
Hessian in such a situation?
- does `optimize.leastsq` do anything different?
- if `x0` is not at a minimum should the covariance matrix be expected to
be positive semi-definite anyway?

_______________________________________________
SciPy-User mailing list
[email protected]
https://mail.python.org/mailman/listinfo/scipy-user