Thanks for your information.
I actually did use second-order method, I thought. If the objective function is the primal function, then I take its gradient, and then used the Jacobian of that gradient.
This is probably true. But the reason might be that I was not using Sparse techniques.
Do you have any idea about this?
I guess I need to invoke some autodiff packages, which works well in conjunction with Sparse techniques. Yet I haven’t figure this out. 
Here is a related link, I’m now reading Advanced tutorial · DifferentiationInterface.jl.