Re: calculating the jacobian for a least-squares problem

> Am 27.03.2018 um 01:57 schrieb Andrew Nelson <[email protected]>:
> 
> 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?
your straightforward approach is ok, especially since you don’t require the highest precision. An alternative would be to use automatic differentiation to calculate the derivatives accurately, e.g. using algopy, theano or tensor flow

> - what is the recommended way (i.e. numerically stable) of inverting the Hessian in such a situation?

If your hess matrix is close to being singular, you could gain some precision by using the QR decomposition of the jacobian. In general to solve a linear system it is recommended to avoid calculating the the inverse matrix.

> - does `optimize.leastsq` do anything different?

leastsq wraps the MINPACK library, which brings it own carefully tuned numeric differentiation routines, and it uses QR decomposition.

> - if `x0` is not at a minimum should the covariance matrix be expected to be positive semi-definite anyway?
If x0 is not a minimum, then there is no guarantee. Even if x0 is a minimum this might by violated due to numerical errors.

best
Gregor

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

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