added optional initial bfgs inverse hessian

This is my first attempt at a pull-request - a really simple additional option for the BFGS method to accept an initial inverse hessian.  I need this feature for a problem in which the standard identity matrix generates errors in the line search routine. 
This allocates a dense ``n x n`` matrix whenever someone calls ``optimize`` for any reason. It should only be allocated when actually needed.