Is it possible to extract the inverse Hessian representation used by L-BFGS after convergence? I would like to use this approximate implicit representation for some further calculations.

I am currently using Optim or NLopt to do the L-BFGS optimization. But I am also open to other packages.

In BFGS yes, see my answer at the bottom here https://stackoverflow.com/questions/43158366/optim-jl-negative-inverse-hessian in LBFGS, no not directly. Why? Because we don’t store it. We calculate the direction based on the history of x and gradient differences. You can access these histories just like in the SO answer and build it yourself. Be aware though, that there can be quite large differences between the hessian approximation and the actual hessian even at “convergence”