ADNLPModels

[PierreMartinon]

@ocots @jbcaillau @gergaud @tmigot @amontoison @gdalle
The issue with Symbolics has been resolved, however the ADNLPModel topic is still relevant, so I propose we continue it here !

[PierreMartinon]

The next step is to use DifferentiationInterface.jl. We could have for a speed-up for performing the AD because I expect the tapes to be better allocated.

But the main speed-up that I expect will be related to the use of the symmetry for the coloring of the Hessian.

If we have less colors, it means less directional derivates. We only use column coloring in ADNLPModels.jl for Hessian and Jacobian right now. I can add it for a release 0.8.2.

This certainly looks promising !
With the old bocop2 I remember that the computation time was roughly divided 50-50 (well, let's say from 1/3-2/3 to 2/3-1/3) between the linear solver part and the NLP problem functions, including derivatives. Among derivatives the Hessian was typically > 90-95% of the total time, apart from a handful of problems with difficult convergence where the Jacobian actually took a significant part.
So I would expect improvements of the Hessian to have the most impact on CTDirect. I'll have to check whether NLPModelIpopt allows acces to the detailed computational times.

[jbcaillau]

@tmigot cute!

[ocots]

@0Yassine0 will try it here: https://github.com/0Yassine0/COTS.jl

[ocots]

@PierreMartinon @jbcaillau how do we handle this? Update of CTDirect.jl or let the user do it?

[PierreMartinon]

Hmm, the typical user may not use these timing statistics, so keep it on a case-by-case basis ?

[jbcaillau]

@ocots @PierreMartinon let's do it the easy way: as we also offer the possibility to retrieve the model (instead of merely solveing), the following is OK (and, yes, in the user's hands):

using OptimalControl
using NLPModelsIpopt
using TimerNLPModels

@def ocp ...
nlp = direct_transcription(ocp; model=:adnlp) # Return an ADNLPModel
nlp = TimerNLPModel(nlp)
stats = ipopt(nlp; print_level = 0)
get_timer(nlp)

[ocots]

Ping @0Yassine0.

[amontoison]

🚀 JuliaSmoothOptimizers/ADNLPModels.jl#258 🚀

[amontoison]

@PierreMartinon @jbcaillau @ocots @0Yassine0 @gergaud
I just released ADNLPModels.jl v0.8.3. 🤘

We use a symmetric coloring now instead of a column coloring for the Hessian of the Lagrangian.
This leads to fewer directional derivatives to compute all nonzeros of the Hessian because we exploit the symmetry.

If you use an optimization solver that needs the Hessian of the Lagrangian (like Ipopt, KNITRO or MadNLP), you will (probably) notice a significant speed-up. In the latest benchmarks, we observed a speed-up greater than x100 🚀 on some problems:
JuliaSmoothOptimizers/ADNLPModels.jl#271 (comment)

These problems likely have a Hessian with an arrowhead form, and exploiting symmetry is particularly important in such cases.
It could highly impact you. 😉

Thanks @gdalle for your implementation of the star coloring 👍
It could be also relevant for you @frapac.

[PierreMartinon]

Hi everyone,
Been working on CTDirect, I did a crude benchmark with 8 problems. Times for increasing grid sizes are below for ADNLPModels 0.8.2 and 0.8.3 (in seconds, excluding precompilation, which takes 50-55s). These runs are for MUMPS with MKL libs. Note that the scaling wrt grid size is not perfect since the number of iterations does change for different grid sizes (eg several problems converged with fewer iterations for N=500 vs N=400, leading to very close overall times).

N	100	200	300	400	500	1000	5000
0.8.2	9	27	115	340	362	?	?
0.8.3	6	9	14	45	46	268	2300
speedup	1.5	3.0	8.2	7.6	7.9	?	?

Not bad at all :D !

[jbcaillau]

thank you @amontoison @gdalle and @PierreMartinon 👍🏽

[gdalle]

With ADNLPModels v0.8.4 (just released), a few improvements to sparsity detection will land right on your doorstep (including reduced compilation time and some array overloads like norm)

[jbcaillau]

well well well. giving it a try 🤞🏾

[amontoison]

@PierreMartinon @jbcaillau @ocots @0Yassine0 @gergaud
I just released ADNLPModels.jl v0.8.7. 🤘 🤘

It provides the opportunity to provide the sparsity pattern of the Jacobians and Hessians
directly if it we already know from the applications or if it was already computed before and dumped in a file.
Tangi did a tutorial about it: https://jso.dev/ADNLPModels.jl/dev/sparsity_pattern/

I also switched to the acycling coloring that I implemented in SparseMatrixColorings.jl for sparse Hessian,
it leads to fewer colors than star or column colorings and thus less directional derivates to recover all non-zeros of the Hessian.
The last benchmark shows a speed-up between 10% and 50%: JuliaSmoothOptimizers/ADNLPModels.jl#294 (comment)

More details about the coloring were also added in the documentation: https://jso.dev/ADNLPModels.jl/dev/sparse/

The last release also uses the new major releases of SparseMatrixColorings.jl (v0.4.0) and SparseConnectivyTracer.jl (v0.6.1)
which lead more speed-ups and robustness.

[jbcaillau]

thanks @amontoison sounds timely and promising 👍🏽🤞🏽!

[ocots]

Thanks for the release and for the tuto!

[jbcaillau]

@PierreMartinon would be nice to re-run your previous tests with the last update
#111 (comment)

[jbcaillau]

JuliaSmoothOptimizers/ADNLPModels.jl#165

[PierreMartinon]

NB. I updated the benchmark and removed 2 rather instable problems.
Runs much faster now, so grid sizes were increased compared to the previous version.

N	250	500	1000	2500	5000	10000
0.8.5	1.6	4.1	11.5	36.2	95.7	271.7
0.8.7	2.0	3.6	7.8	29.3	79.1	197.2
speedup	0.8	1.14	1.47	1.24	1.21	1.38

Seems we got a 20-40% speed gain, nice !

[jbcaillau]

good, thanks @PierreMartinon 👍🏽

[amontoison]

New benchmarks on the ACOPF problems:

• https://discourse.julialang.org/t/ac-optimal-power-flow-in-various-nonlinear-optimization-frameworks/78486/98

I think that replacing ReverseDiff.jl with Enzyme.jl will be the easiest way to gain a speed-up in the future.
We have only updated sparsity detection and coloring since last year.

[gdalle]

I'm currently working on multiple arguments for DifferentiationInterface, which was the last hurdle for the hessian of the lagrangian with a fixed multiplicator. Once that lands, I think we can consider the switch to DI

[amontoison]

I will discuss about the Hessian of the Lagrangian with William Moses during AD2024.
Maybe we can optimize things at low-level.

[gdalle]

In that case such optimizations are likely to land either in Enzyme or in DI's Enzyme extension 😉

[PierreMartinon]

Allright, benchmark with the new CTModels (abstract) !

Trapeze

julia> bench()
Grid size    250: time (s) =    0.5
Grid size    500: time (s) =    1.4
Grid size   1000: time (s) =    3.5
Grid size   2500: time (s) =   17.3
Grid size   5000: time (s) =   71.0

julia> bench(;adnlp_backend=:manual)
Grid size    250: time (s) =    0.8
Grid size    500: time (s) =    1.7
Grid size   1000: time (s) =    3.5
Grid size   2500: time (s) =   11.3
Grid size   5000: time (s) =   24.0

Midpoint

julia> bench(;disc_method=:midpoint)
Grid size    250: time (s) =    0.7
Grid size    500: time (s) =    1.6
Grid size   1000: time (s) =    4.1
Grid size   2500: time (s) =   17.9
Grid size   5000: time (s) =   62.2

julia> bench(;disc_method=:midpoint, adnlp_backend=:manual)
Grid size    250: time (s) =    0.8
Grid size    500: time (s) =    1.9
Grid size   1000: time (s) =    3.8
Grid size   2500: time (s) =   11.6
Grid size   5000: time (s) =   24.2

Gauss Legendre 2

julia> bench(;disc_method=:gauss_legendre_2)
Grid size    250: time (s) =    2.7
Grid size    500: time (s) =    9.1
Grid size   1000: time (s) =   22.5
<grid size 2500 fails one of the tests -_->
Grid size   5000: time (s) =  366.6

julia> bench(;disc_method=:gauss_legendre_2, adnlp_backend=:manual)
Grid size    250: time (s) =    3.3
Grid size    500: time (s) =    8.0
Grid size   1000: time (s) =   18.5
Grid size   2500: time (s) =   55.5
Grid size   5000: time (s) =  131.0

Gauss Legendre 3

julia> bench(;disc_method=:gauss_legendre_3)
Grid size    250: time (s) =    5.1
Grid size    500: time (s) =   15.9
Grid size   1000: time (s) =   34.4
Grid size   2500: time (s) =  160.4
Grid size   5000: time (s) =  624.0

julia> bench(;disc_method=:gauss_legendre_3, adnlp_backend=:manual)
Grid size    250: time (s) =    6.3
Grid size    500: time (s) =   15.8
Grid size   1000: time (s) =   38.7
Grid size   2500: time (s) =  104.9
Grid size   5000: time (s) =  194.1

[jbcaillau]

super nice @PierreMartinon 👍🏽 extra wish 🎅🏾:

• would like to see allocations 🙏🏽
• see that trapeze is always slightly faster than its dual midpoint scheme (with the advantage that trapeze has the same gridpoints)
• for the same optim tolerance, higher order schemes are slower but is there a gain at least in the number of iterations? (= would be nice to display time + alloc + iteration number 🥲)

[jbcaillau]

@rand-asswad please give a try to the new high order schemes on your problem 🤞🏽