Issue recovering inverse Hessian from minimize()

[rickecon]

The inverse Hessian is a nice estimator for the variance-covariance matrix of MLE parameter estimates (see MLE notebook discussion in Section 4). There are two problems with this.

1. Only one of the three constrained minimizer methods ('L-BFGS-B') return the inverse Hessian.
2. Sometimes, the 'L-BFGS-B' method returns the inverse Hessian as a <KxK LbfgsInvHessProduct with dtype=float64>, where K is the number of parameters being estimated.

This issue presents some solutions to these problems as well as some unresolved questions. @merubhanot @dpzhang @olivianatan @alejandroparraguez @bobaekang @ykim17 @kkost84

[rickecon]

With regard to problem 1, suppose that it is either the 'TNC' or 'SLSQP' constrained minimizer method that successfully converges to an answer. A cheap work around is to use the estimated parameter values from 'TNC' or 'SLSQP' as initial values for 'L-BFGS-B' and then recover the estimate of inverse Hessian from the 'L-BFGS-B' estimates. Make sure that the parameter estimates and the log likelihood function don't change much in this second step. You can also sometimes get the unconstrained minimizer to converge by using as starting values the estimates from either 'TNC' or 'SLSQP'. The unconstrained minimizer also returns the estimate of the inverse Hessian.

[rickecon]

With regard to problem 2, I have not yet figured out how to recover the hess_inv object when it is stored as a<KxK LbfgsInvHessProduct with dtype=float64>, where K is the number of parameters being estimated. However, I looked in the minimize() documentation in the section on what the function returns. This led me to the OptimizeResult documentation. The description of the hess_inv object under attributes says, "The type of this attribute may be either np.ndarray or scipy.sparse.linalg.LinearOperator."

The scipy.sparse.linalg.LinearOperator documentation is complicated, but I think the answer might be to use the .matvec() method. For example, suppose one performed a constrained minimization in which 5 parameters were chosen to minimize the criterion function such that the hess_inv object were stored as a <5x5 LbfgsInvHessProduct with dtype=float64> object in the OptimizeResult object named results. I think the inverse Hessian could be recovered by typing the following code.

VCV = results.hess_inv.matmat(np.eye(5))

In other words, you perform a matrix multiplication of the hess_inv object by the identity matrix of the same size to return the original matrix. I am not 100% sure this is right, nor do I understand why this would be the case. But I think it might be correct.

[merubhanot]

I ended up using the following .todense() command to obtain the inverse hessian. Example code:

results.hess_inv.todense()

[kkost84]

I also used a different method that is somewhat unclear in how it works but after checking it returns the same thing as Professor Evans matmat method and also Meru's todense method.

results.hess_inv*np.identity(K) where K is the dimension of the operator