As the title says, I would like some insight on what Julia Optimization package would be most appropriate to efficiently solve nonlinear problems with potentially exploitable sparsity structures (e.g. block diagonal Hessian).

I have come across various packages:

NLPModels.jl emphasises the capability for user-defined jac/hess structure;

JuMP.jl, Optim.jl, NLOpt.jl all seem to support nonlinear problems, but I’m not sure what happens under the hood for exploiting sparsity.
Perhaps someone with experience on this could help?

It relies on solvers that do not build Hessian directly (it uses for example L-BFGS) but it can handle problems with sparse cost functions such as LASSO.

For user-defined functions, we disable all hessians, and compute the gradient with forward mode AD. Depending on the model, this can cause significant issues.

If at all possible, you should write out your constraints algebraically.