Re: calculating the jacobian for a least-squares problem

On Mon, Mar 26, 2018 at 7:57 PM, Andrew Nelson <[email protected]> wrote:
> 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 never looked much at the details of MCMC.
But if your data or starting point doesn't provide good information about the
Hessian, then, I think, you could shrink the hessian to or combine it with the
prior covariance matrix, e.g. use a weighted average.

Josef


 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
>