Re: calculating the jacobian for a least-squares problem

On Wed, Mar 28, 2018 at 7:33 PM, Andrew Nelson <[email protected]> wrote:
> I'm using the Hessian to calculate the covariance matrix for parameter
> estimates in least squares, i.e. the equivalent of `pcov` in `curve_fit` (I
> don't want to do a fit, I just want the covariance around the current
> location).

But if you are just using it to initialize a Bayesian estimator, then you
don't need the covariance from the data by itself.

What I guess:
You are using the cross-product of the jacobian (outer product of gradient).

This should in general be positive definite.
If it is not positive definite, then locally at least one of the
parameters is not
identified, i.e. your X in analogy to linear regression is singular.
You can fiddle with numerical precision to maybe get it noisily PSD, and
then you have just an almost singular jacobian with a very noise inverse.

pinv(jac)

Josef

>
> On 29 March 2018 at 03:05, <[email protected]> wrote:
>>
>> 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
>> >
>> _______________________________________________
>> SciPy-User mailing list
>> [email protected]
>> https://mail.python.org/mailman/listinfo/scipy-user
>
>
>
>
> --
> _____________________________________
> Dr. Andrew Nelson
>
>
> _____________________________________
>
> _______________________________________________
> SciPy-User mailing list
> [email protected]
> https://mail.python.org/mailman/listinfo/scipy-user
>