LBFGS for use in large maximum likelihood problem

I've got another set of questions about LBFGS, maximum likelihood, and 
numerical issues.

I'm using the inverse of the hessian of the log-likelihood function 
evaluated at the maximum likelihood estimates (from LBFGS) to 
approximate the uncertainties in the estimated coefficients.

Calculating the hessian directly from the log-likelihood expression 
takes thousands of seconds... If i calculate the hessian of one symbolic 
term of the sum, apply the "optimize" function, and then numerically 
evaluate the sum of the optimized expression it takes a few seconds (8 
or so). However optimize says

"`optimize' has met up with a special form - answer may be wrong."

I checked the result against a non-optimized version (which takes 500 
seconds or so) and I get the same numerical values.

When should I worry about relying on the "optimize" results?

Also, can I use optimize on the figure of merit expression for LBFGS to 
get faster results? Is there a way to hand LBFGS an expression for the 
gradient (also optimized)? The manual doesn't mention anything, but I 
thought there was some work done on this?

Thanks much,
Dan