Multi-precision and LBFGS

[edljk]

Following the multi-precision example of the documentation I tried a quasi-Newton call without success:

using MadNLP, ExaModels, Quadmath

function airport_model(T)
    N = 42
    # Data
    r = T[0.09 , 0.3, 0.09, 0.45, 0.5, 0.04, 0.1, 0.02, 0.02, 0.07, 0.4, 0.045, 0.05, 0.056, 0.36, 0.08, 0.07, 0.36, 0.67, 0.38, 0.37, 0.05, 0.4, 0.66, 0.05, 0.07, 0.08, 0.3, 0.31, 0.49, 0.09, 0.46, 0.12, 0.07, 0.07, 0.09, 0.05, 0.13, 0.16, 0.46, 0.25, 0.1]
    cx = T[-6.3, -7.8, -9.0, -7.2, -5.7, -1.9, -3.5, -0.5, 1.4, 4.0, 2.1, 5.5, 5.7, 5.7, 3.8, 5.3, 4.7, 3.3, 0.0, -1.0, -0.4, 4.2, 3.2, 1.7, 3.3, 2.0, 0.7, 0.1, -0.1, -3.5, -4.0, -2.7, -0.5, -2.9, -1.2, -0.4, -0.1, -1.0, -1.7, -2.1, -1.8, 0.0]
    cy = T[8.0, 5.1, 2.0, 2.6, 5.5, 7.1, 5.9, 6.6, 6.1, 5.6, 4.9, 4.7, 4.3, 3.6, 4.1, 3.0, 2.4, 3.0, 4.7, 3.4, 2.3, 1.5, 0.5, -1.7, -2.0, -3.1, -3.5, -2.4, -1.3, 0.0, -1.7, -2.1, -0.4, -2.9, -3.4, -4.3, -5.2, -6.5, -7.5, -6.4, -5.1, 0.0]
    # Wrap all data in a single iterator for ExaModels
    data = [(i, cx[i], cy[i], r[i]) for i in 1:N]
    IJ = [(i, j) for i in 1:N-1 for j in i+1:N]
    # Write model using ExaModels
    core = ExaModels.ExaCore(T)
    x = ExaModels.variable(core, 1:N, lvar = -10.0, uvar=10.0)
    y = ExaModels.variable(core, 1:N, lvar = -10.0, uvar=10.0)
    ExaModels.objective(
        core,
        ((x[i] - x[j])^2 + (y[i] - y[j])^2) for (i, j) in IJ
    )
    ExaModels.constraint(core, (x[i]-dcx)^2 + (y[i] - dcy)^2 - dr for (i, dcx, dcy, dr) in data; lcon=-Inf)
    return ExaModels.ExaModel(core)
end

nlp_128 = airport_model(Float128)
results_128 = madnlp(nlp_128; linear_solver = LDLSolver,
                                        hessian_approximation = MadNLP.CompactLBFGS)

ERROR: MethodError: no method matching trsm!(::Char, ::Char, ::Char, ::Char, ::Float128, ::Matrix{Float128}, ::SubArray{Float128, 2, Matrix{…}, Tuple{…}, true})
The function `trsm!` exists, but no method is defined for this combination of argument types.
...

[amontoison]

@edljk LBFGS is under a major revamp.
I optimized it to use BLAS and LAPACK as much as possible (single / double support) but I need to add some generic fallback of BLAS / LAPACK foutines for non-standard precision.

I will try to fix it for version 0.9 of MadNLP.
Note that it should work in Float32.

[edljk]

Perfect. Thank you for the great package!

Btw with Float32 I get another difficulty

 1  1.3630479e+04 2.59e+01 3.65e+02 3.69e+00  -1.0 4.50e+00    -  6.87e-01 1.00e+00h  1
   2r 1.3630479e+04 9.54e-07 9.99e+02 2.59e+01   1.4 2.35e-10    -  6.87e-01 1.00e+00h  1
EXIT: Restoration Failed
"Execution stats: Restoration Failed"

[amontoison]

Ok, good to know.
I will cross-reference this issue with #537.
Did you try it with the branch master of MadNLP?

[edljk]

Yes the default branch

(@v1.12) pkg> status MadNLP 
Status `~/.julia/environments/v1.12/Project.toml`
  [2621e9c9] MadNLP v0.8.12

[amontoison]

@edljk I wanted to say the current master branch with the modifications of the last 4 months that we didn't released.

pkg> add MadNLP#master

[edljk]

Oups, sorry.
Switching to master, Float32 and MadNLP.CompactLBFGS works fine up to tol = 1e-3 using LapackCPUSolver. On the other hand LDLSolver (perhaps not relevant in this context..) an error:

0  0.0000000e+00 1.04e+02 9.98e-01 1.00e+01  -1.0     -   0.00e+00  1  0 
   1  1.3630478e+04 2.59e+01 3.65e+02 3.69e+00  -1.0     -   1.00e+00  1  1h
EXIT: Internal Error.
ERROR: MethodError: no method matching solve!(::MadNLP.MumpsSolver{Float32}, ::Vector{Float64})
The function `solve!` exists, but no method is defined for this combination of argument types.

[amontoison]

We should have Float32 everywhere. It is bug in MadNLP 🙂
I will fix it for the release 0.9.

[amontoison]

@edljk I fixed the issue for single precision in #567.
For quadruple precision, I need to add a generic fallback for these functions:
https://github.com/MadNLP/MadNLP.jl/blob/master/src/utils.jl#L40-L65

Do you need the support for Float128 for your problem or it was just to test LBFGS?
It will help to decide what is the priority of this task.

[amontoison]

using MadNLP, ExaModels, Quadmath

function airport_model(T)
    N = 42
    # Data
    r = T[0.09 , 0.3, 0.09, 0.45, 0.5, 0.04, 0.1, 0.02, 0.02, 0.07, 0.4, 0.045, 0.05, 0.056, 0.36, 0.08, 0.07, 0.36, 0.67, 0.38, 0.37, 0.05, 0.4, 0.66, 0.05, 0.07, 0.08, 0.3, 0.31, 0.49, 0.09, 0.46, 0.12, 0.07, 0.07, 0.09, 0.05, 0.13, 0.16, 0.46, 0.25, 0.1]
    cx = T[-6.3, -7.8, -9.0, -7.2, -5.7, -1.9, -3.5, -0.5, 1.4, 4.0, 2.1, 5.5, 5.7, 5.7, 3.8, 5.3, 4.7, 3.3, 0.0, -1.0, -0.4, 4.2, 3.2, 1.7, 3.3, 2.0, 0.7, 0.1, -0.1, -3.5, -4.0, -2.7, -0.5, -2.9, -1.2, -0.4, -0.1, -1.0, -1.7, -2.1, -1.8, 0.0]
    cy = T[8.0, 5.1, 2.0, 2.6, 5.5, 7.1, 5.9, 6.6, 6.1, 5.6, 4.9, 4.7, 4.3, 3.6, 4.1, 3.0, 2.4, 3.0, 4.7, 3.4, 2.3, 1.5, 0.5, -1.7, -2.0, -3.1, -3.5, -2.4, -1.3, 0.0, -1.7, -2.1, -0.4, -2.9, -3.4, -4.3, -5.2, -6.5, -7.5, -6.4, -5.1, 0.0]
    # Wrap all data in a single iterator for ExaModels
    data = [(i, cx[i], cy[i], r[i]) for i in 1:N]
    IJ = [(i, j) for i in 1:N-1 for j in i+1:N]
    # Write model using ExaModels
    core = ExaModels.ExaCore(T)
    x = ExaModels.variable(core, 1:N, lvar = -10.0, uvar=10.0)
    y = ExaModels.variable(core, 1:N, lvar = -10.0, uvar=10.0)
    ExaModels.objective(
        core,
        ((x[i] - x[j])^2 + (y[i] - y[j])^2) for (i, j) in IJ
    )
    ExaModels.constraint(core, (x[i]-dcx)^2 + (y[i] - dcy)^2 - dr for (i, dcx, dcy, dr) in data; lcon=-Inf)
    return ExaModels.ExaModel(core)
end

nlp_32 = airport_model(Float32)
results_32 = madnlp(nlp_32; hessian_approximation = MadNLP.CompactLBFGS)

julia> results_32 = madnlp(nlp_32; hessian_approximation = MadNLP.CompactLBFGS)
This is MadNLP version v0.8.12, running with MUMPS v5.8.2

Number of nonzeros in constraint Jacobian............:       84
Number of nonzeros in Lagrangian Hessian.............:     5250

Total number of variables............................:       84
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       84
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:       42
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       42

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  0.0000000e+00 1.04e+02 9.98e-01 1.00e+01  -1.0     -   0.00e+00  1  0 
   1  1.3630478e+04 2.59e+01 3.65e+02 3.69e+00  -1.0     -   1.00e+00  1  1h
   2  1.3671312e+04 2.56e+01 3.60e+02 1.35e+00  -1.0     -   1.26e-02  1  1h
   3  1.5200467e+04 2.25e+01 3.15e+02 2.68e+00  -1.0     -   1.25e-01  1  1h
   4  1.7169811e+04 1.94e+01 2.71e+02 2.48e+00  -1.0     -   1.45e-01  1  1h
   5  1.8580420e+04 1.74e+01 2.45e+02 2.47e+00  -1.0     -   1.06e-01  1  1h
   6  3.0672289e+04 5.99e+00 2.09e+02 1.53e+00  -1.0     -   8.24e-01  1  1h
   7  3.6431402e+04 2.97e+00 1.67e+02 6.21e-01  -1.0     -   5.89e-01  1  1h
   8  4.3683820e+04 7.23e-01 2.08e+02 1.44e-01  -1.0     -   1.00e+00  1  1h
   9  4.5871480e+04 2.85e-01 2.16e+02 1.06e-01  -1.0     -   7.86e-01  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  10  4.5908656e+04 2.72e-01 1.91e+02 1.29e-01  -1.0     -   4.54e-02  1  1h
  11  4.6929359e+04 1.15e-01 1.04e+02 2.17e-01  -1.0     -   1.00e+00  1  1h
  12  4.7714723e+04 3.05e-02 7.19e+01 9.37e-02  -1.0     -   1.00e+00  1  1h
  13  4.7915961e+04 5.62e-03 3.17e+01 8.15e-02  -1.0     -   1.00e+00  1  1h
  14  4.7951871e+04 1.05e-03 1.13e+01 7.86e-02  -1.0     -   1.00e+00  1  1h
  15  4.7956941e+04 6.53e-05 2.57e+00 7.78e-02  -1.0     -   1.00e+00  1  1h
  16  4.7957004e+04 9.07e-06 2.15e+00 7.77e-02  -1.0     -   1.00e+00  1  1h
  17  4.7956730e+04 2.72e-04 8.26e+00 7.76e-02  -1.0     -   1.00e+00  1  1h
  18  4.7956719e+04 8.40e-05 6.15e+00 7.77e-02  -1.0     -   1.00e+00  1  1h
  19  4.7956848e+04 2.42e-05 2.84e-01 7.76e-02  -1.0     -   1.00e+00  1  1h
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  20  4.7953512e+04 2.23e-06 6.02e-01 1.73e-02  -1.7     -   1.00e+00  1  1f
  21  4.7953484e+04 4.81e-07 3.95e-01 1.55e-02  -1.7     -   1.00e+00  1  1h
  22  4.7953484e+04 2.63e-07 5.66e-02 1.55e-02  -1.7     -   1.00e+00  1  1h
  23  4.7952848e+04 1.99e-07 3.90e-02 2.23e-03  -2.5     -   1.00e+00  2  1h
  24  4.7952832e+04 2.96e-07 4.96e-02 2.20e-03  -2.5     -   1.00e+00  1  1h
  25  4.7952824e+04 2.92e-07 2.74e-02 2.20e-03  -2.5     -   2.50e-01  1  3h
  26  4.7952711e+04 3.05e-07 1.31e-02 1.17e-04  -3.8     -   1.00e+00  2  1h
  27  4.7952727e+04 2.38e-07 8.33e-03 1.17e-04  -3.8     -   5.00e-01  1  2h
  28  4.7952727e+04 2.98e-07 8.20e-03 1.17e-04  -3.8     -   3.12e-02  1  6h
  29  4.7952727e+04 2.98e-07 7.14e+02 1.87e+01  -3.8     -   1.56e-02  1  7 
iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
  30  4.7952730e+04 2.89e-07 2.17e+00 4.12e-02  -3.8     -   1.00e+00  1  1h
  31  4.7952727e+04 2.38e-07 6.64e-03 1.17e-04  -3.8     -   1.00e+00  1  1h
  32  4.7952723e+04 4.22e-07 2.25e-03 1.17e-04  -3.8     -   1.00e+00  1  1h
  33  4.7952723e+04 2.44e-07 2.12e-03 1.17e-04  -3.8     -   6.25e-02  1  5h
  34  4.7952723e+04 2.44e-07 7.14e+02 1.87e+01  -3.8     -   3.12e-02  1  6 
  35  4.7952727e+04 2.83e-07 2.17e+00 4.12e-02  -3.8     -   1.00e+00  1  1h
  36  4.7952730e+04 2.24e-07 1.39e-03 1.17e-04  -3.8     -   1.00e+00  1  1h
  37  4.7952719e+04 1.66e-07 2.21e-03 7.78e-06  -5.0     -   1.00e+00  2  1h
  38  4.7952719e+04 4.23e-07 1.07e-03 7.77e-06  -5.0     -   1.00e+00  1  1H
  39  4.7952723e+04 1.51e-07 9.60e-04 7.77e-06  -5.0     -   1.00e+00  1  1H

Number of Iterations....: 39

                                   (scaled)                 (unscaled)
Objective...............:   4.7952722656250000e+04    4.7952722656250000e+04
Dual infeasibility......:   9.6011284040287137e-04    9.6011284040287137e-04
Constraint violation....:   1.5069905145992379e-07    1.5069905145992379e-07
Complementarity.........:   7.7717422755085863e-06    7.7717422755085863e-06
Overall NLP error.......:   9.6011284040287137e-04    9.6011284040287137e-04

Number of objective function evaluations              = 67
Number of objective gradient evaluations              = 38
Number of constraint evaluations                      = 67
Number of constraint Jacobian evaluations             = 40
Number of Lagrangian Hessian evaluations              = 0

Total wall secs in initialization                     =  0.001
Total wall secs in linear solver                      =  0.035
Total wall secs in NLP function evaluations           =  0.001
Total wall secs in solver (w/o init./fun./lin. alg.)  =  2.225
Total wall secs                                       =  2.262

EXIT: Optimal Solution Found (tol = 1.0e-03).
"Execution stats: Optimal Solution Found (tol = 1.0e-03)."

[edljk]

Thank you for the quick response!!
That's very strange, I just updated MadNLP, copied and pasted the Float32 example and got

julia> results_32 = madnlp(nlp_32; hessian_approximation = MadNLP.CompactLBFGS)
This is MadNLP version v0.8.12, running with MUMPS v5.8.2

Number of nonzeros in constraint Jacobian............:       84
Number of nonzeros in Lagrangian Hessian.............:     5250

Total number of variables............................:       84
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       84
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:       42
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:       42

iter    objective    inf_pr   inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
   0  0.0000000e+00 1.04e+02 9.98e-01 1.00e+01  -1.0     -   0.00e+00  1  0 
   1  1.3630478e+04 2.59e+01 3.65e+02 3.69e+00  -1.0     -   1.00e+00  1  1h
EXIT: Internal Error.
ERROR: MethodError: no method matching solve!(::MumpsSolver{Float32}, ::Vector{Float64})
The function `solve!` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  solve!(::LapackCPUSolver{Float64}, ::Vector{Float64})
   @ MadNLP ~/.julia/packages/MadNLP/X57EO/src/LinearSolvers/lapack.jl:82
  solve!(::LapackCPUSolver, ::AbstractVector)

[amontoison]

@edljk You need the branch am/float32_lbfgs.
My fix will be merged tomorrow by François (@frapac)
The related PR is #567.

[edljk]

Great!
Yes, I am very interested in the multi precision support (Float128, Float64x2, Float64x4..) for quasi Newton methods in NLP. I still have a few steps to complete in the coming days before moving to optimization, so no rush ;-)

[edljk]

I'm back working on multi-precision optimization. In cases where an analytic Hessian is provided, what types are supported? Float64x2 and Float64x4 don't appear to be supported by the linear solvers.

[sshin23]

Below is another option:

madnlp(m; kkt_system=MadNLP.SparseCondensedKKTSystem,linear_solver=MadNLP.LDLSolver, tol=1e-16)

LDLFactorization.jl is a pure julia solver, so it will work with whatever primitive fp type you may have. kkt_system=MadNLP.SparseCondensedKKTSystem is necessary to make numerical linear algebra compatible with LDLFactorization.