Running MadNLP in arbitrary precision
MadNLP supports solving optimization problems in arbitrary precision. By default, MadNLP adapts its precision following the NLPModel passed in input. Most models are using Float64 (in fact, almost all optimization modelers are implemented using double precision), but for certain applications it can be useful to use arbitrary precision to get more accurate solution.
There exists different packages to instantiate a optimization model in arbitrary precision in Julia. Most of them leverage the flexibility offered by NLPModels.jl. In particular, we recommend:
Defining a problem in arbitrary precision
As a demonstration, we implement the model airport from CUTEst using ExaModels. The code writes:
using ExaModels
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)
endairport_model (generic function with 1 method)The function airport_model takes as input the type used to define the model in ExaModels. For example, ExaCore(Float64) instantiates a model with Float64, whereas ExaCore(Float32) instantiates a model using Float32. Thus, instantiating the instance airport using Float32 simply amounts to
nlp = airport_model(Float32)An ExaModel{Float32, Vector{Float32}, ...}
Problem name: Generic
All variables: ████████████████████ 84 All constraints: ████████████████████ 42
free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ████████████████████ 42
low/upp: ████████████████████ 84 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: (-47.06% sparsity) 5250 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ████████████████████ 42
nnzj: ( 97.62% sparsity) 84
lin_nnzj: (------% sparsity)
nln_nnzj: ( 97.62% sparsity) 84
We verify that the model is correctly instantiated using Float32:
x0 = NLPModels.get_x0(nlp)
println(typeof(x0))Vector{Float32}Solving a problem in Float32
Now that we have instantiated our model in Float32, we solve it using MadNLP. As nlp is using Float32, MadNLP will automatically adjust its internal types to Float32 during the instantiation. By default, the convergence tolerance is also adjusted to the input type, such that tol = sqrt(eps(T)). Hence, in our case the tolerance is set automatically to
tol = sqrt(eps(Float32))0.00034526698f0We solve the problem using Lapack as linear solver:
results = madnlp(nlp; linear_solver=LapackCPUSolver)
nothingThis is MadNLP version v0.9.1, running with Lapack-CPU (BUNCHKAUFMAN)
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.4805764e+01 1.00e+02 2.63e+01 6.04e+00 -1.0 - 3.37e-02 1 1h
2 1.3660695e+04 2.50e+01 1.79e+02 1.87e+00 -1.0 - 1.00e+00 1 1h
3 1.5777956e+04 2.10e+01 1.55e+02 8.18e-01 -1.0 - 1.68e-01 1 1h
4 2.2549984e+04 1.21e+01 1.30e+02 3.48e-01 -1.0 - 4.83e-01 1 1h
5 3.6154164e+04 3.00e+00 2.13e+02 2.72e-01 -1.0 - 1.00e+00 1 1h
6 4.0256086e+04 1.54e+00 1.64e+02 1.19e-01 -1.0 - 5.66e-01 1 1h
7 4.5447668e+04 3.63e-01 1.99e+02 1.15e-01 -1.0 - 1.00e+00 1 1h
8 4.7406473e+04 7.26e-02 1.00e+02 1.00e-01 -1.0 - 1.00e+00 1 1h
9 4.7896695e+04 8.12e-03 3.88e+01 9.06e-02 -1.0 - 1.00e+00 1 1h
iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
10 4.7955184e+04 3.16e-04 6.32e+00 7.81e-02 -1.0 - 1.00e+00 1 1h
11 4.7956895e+04 1.32e-06 6.31e-02 7.79e-02 -1.0 - 1.00e+00 1 1h
12 4.7952863e+04 3.40e-06 1.26e-03 2.85e-03 -2.5 - 1.00e+00 1 1f
13 4.7952723e+04 3.04e-07 2.99e-04 8.49e-06 -5.0 - 1.00e+00 2 1h
Number of Iterations....: 13
(scaled) (unscaled)
Objective...............: 4.7952722656250000e+04 4.7952722656250000e+04
Dual infeasibility......: 2.9884013929404318e-04 2.9884013929404318e-04
Constraint violation....: 3.0411138141062111e-07 3.0411138141062111e-07
Complementarity.........: 8.4870234786649235e-06 8.4870234786649235e-06
Overall NLP error.......: 2.9884013929404318e-04 2.9884013929404318e-04
Number of objective function evaluations = 14
Number of objective gradient evaluations = 14
Number of constraint evaluations = 14
Number of constraint Jacobian evaluations = 14
Number of Lagrangian Hessian evaluations = 13
Number of KKT factorizations = 14
Number of KKT backsolves = 15
Total wall secs in initialization = 4.198 s
Total wall secs in linear solver = 0.004 s
Total wall secs in NLP function evaluations = 0.001 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 20.398 s
Total wall secs = 24.602 s
EXIT: Optimal Solution Found (tol = 1.0e-03).Note that the distribution of Lapack shipped with Julia supports Float32, so here we do not have to worry whether the type is supported by the linear solver. Almost all linear solvers shipped with MadNLP supports Float32.
The final solution is
results.solution84-element Vector{Float32}:
-6.1042237
-7.3201675
-8.7022505
-6.547783
-5.1567717
-1.8469505
-3.3253756
-0.49251527
1.3583232
3.8234024
⋮
-2.6840556
-3.1422696
-4.000223
-4.9764223
-6.141531
-7.106216
-5.7435718
-4.6168227
0.30903557and the objective is
results.objective47952.723f0For completeness, we compare with the solution returned when we solve the same problem using Float64:
nlp_64 = airport_model(Float64)
results_64 = madnlp(nlp_64; linear_solver=LapackCPUSolver)
nothingThis is MadNLP version v0.9.1, running with Lapack-CPU (BUNCHKAUFMAN)
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.4805731e+01 1.00e+02 2.63e+01 6.04e+00 -1.0 - 3.37e-02 1 1h
2 1.3660700e+04 2.50e+01 1.79e+02 1.87e+00 -1.0 - 1.00e+00 1 1h
3 1.5777961e+04 2.10e+01 1.55e+02 8.18e-01 -1.0 - 1.68e-01 1 1h
4 2.2549946e+04 1.21e+01 1.30e+02 3.48e-01 -1.0 - 4.83e-01 1 1h
5 3.6154159e+04 3.00e+00 2.13e+02 2.72e-01 -1.0 - 1.00e+00 1 1h
6 4.0256063e+04 1.54e+00 1.64e+02 1.19e-01 -1.0 - 5.66e-01 1 1h
7 4.5447681e+04 3.63e-01 1.99e+02 1.15e-01 -1.0 - 1.00e+00 1 1h
8 4.7406467e+04 7.26e-02 1.00e+02 1.00e-01 -1.0 - 1.00e+00 1 1h
9 4.7896713e+04 8.12e-03 3.88e+01 9.06e-02 -1.0 - 1.00e+00 1 1h
iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
10 4.7955178e+04 3.16e-04 6.32e+00 7.81e-02 -1.0 - 1.00e+00 1 1h
11 4.7956896e+04 1.29e-06 6.32e-02 7.79e-02 -1.0 - 1.00e+00 1 1h
12 4.7952823e+04 3.39e-06 1.21e-03 2.85e-03 -2.5 - 1.00e+00 1 1f
13 4.7952701e+04 4.21e-09 1.24e-06 2.16e-06 -5.7 - 1.00e+00 1 1h
14 4.7952701e+04 2.18e-15 7.42e-13 1.95e-09 -8.6 - 1.00e+00 2 1h
Number of Iterations....: 14
(scaled) (unscaled)
Objective...............: 4.7952701409727124e+04 4.7952701409727124e+04
Dual infeasibility......: 7.4217751281706175e-13 7.4217751281706175e-13
Constraint violation....: 2.1754946989537970e-15 2.1754946989537970e-15
Complementarity.........: 1.9478836876136048e-09 1.9478836876136048e-09
Overall NLP error.......: 1.9478836876136048e-09 1.9478836876136048e-09
Number of objective function evaluations = 15
Number of objective gradient evaluations = 15
Number of constraint evaluations = 15
Number of constraint Jacobian evaluations = 15
Number of Lagrangian Hessian evaluations = 14
Number of KKT factorizations = 15
Number of KKT backsolves = 16
Total wall secs in initialization = 2.757 s
Total wall secs in linear solver = 0.006 s
Total wall secs in NLP function evaluations = 0.001 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 8.182 s
Total wall secs = 10.946 s
EXIT: Optimal Solution Found (tol = 1.0e-08).The final objective is now
results_64.objective47952.701409727124As expected when solving an optimization problem with Float32, the relative difference between the two solutions is far from being negligible:
rel_diff = abs(results.objective - results_64.objective) / results_64.objective4.4307249125524247e-7Solving a problem in Float128
Now, we go in the opposite direction and solve a problem using Float128 to get a better accuracy. We start by loading the library Quadmath to work with quadruple precision:
using QuadmathWe can instantiate our problem using Float128 directly as:
nlp_128 = airport_model(Float128)An ExaModel{Quadmath.Float128, Vector{Quadmath.Float128}, ...}
Problem name: Generic
All variables: ████████████████████ 84 All constraints: ████████████████████ 42
free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 free: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 lower: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
upper: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 upper: ████████████████████ 42
low/upp: ████████████████████ 84 low/upp: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 fixed: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0 infeas: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nnzh: (-47.06% sparsity) 5250 linear: ⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅ 0
nonlinear: ████████████████████ 42
nnzj: ( 97.62% sparsity) 84
lin_nnzj: (------% sparsity)
nln_nnzj: ( 97.62% sparsity) 84
On the contrary to Float32, a few linear solvers support Float128 out of the box. Currently, the only solvers suporting quadruple in MadNLP are LDLSolver and the HSL solvers (require MadNLPHSL). LDLSolvers uses an LDL factorization implemented in pure Julia. The solver LDLSolver is not adapted to solve large-scale nonconvex nonlinear programs, but works if the problem is small enough (as it is the case here).
Once we have replaced the solver by LDLSolver, solving the problem with MadNLP in Float128 just amounts to
results_128 = madnlp(nlp_128; linear_solver=LDLSolver)
nothingThis is MadNLP version v0.9.1, running with LDLFactorizations v0.10.1
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.4805731e+01 1.00e+02 2.63e+01 6.04e+00 -1.0 - 3.37e-02 1 1h
2 1.3660700e+04 2.50e+01 1.79e+02 1.87e+00 -1.0 - 1.00e+00 1 1h
3 1.5777961e+04 2.10e+01 1.55e+02 8.18e-01 -1.0 - 1.68e-01 1 1h
4 2.2549946e+04 1.21e+01 1.30e+02 3.48e-01 -1.0 - 4.83e-01 1 1h
5 3.6154159e+04 3.00e+00 2.13e+02 2.72e-01 -1.0 - 1.00e+00 1 1h
6 4.0256063e+04 1.54e+00 1.64e+02 1.19e-01 -1.0 - 5.66e-01 1 1h
7 4.5447681e+04 3.63e-01 1.99e+02 1.15e-01 -1.0 - 1.00e+00 1 1h
8 4.7406467e+04 7.26e-02 1.00e+02 1.00e-01 -1.0 - 1.00e+00 1 1h
9 4.7896713e+04 8.12e-03 3.88e+01 9.06e-02 -1.0 - 1.00e+00 1 1h
iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls
10 4.7955178e+04 3.16e-04 6.32e+00 7.81e-02 -1.0 - 1.00e+00 1 1h
11 4.7956896e+04 1.29e-06 6.32e-02 7.79e-02 -1.0 - 1.00e+00 1 1h
12 4.7952823e+04 3.39e-06 1.21e-03 2.85e-03 -2.5 - 1.00e+00 1 1f
13 4.7952701e+04 4.21e-09 1.24e-06 2.16e-06 -5.7 - 1.00e+00 1 1h
14 4.7952701e+04 2.17e-15 5.84e-13 1.95e-09 -8.6 - 1.00e+00 1 1h
15 4.7952701e+04 2.31e-21 8.16e-19 1.15e-18 -18.0 - 1.00e+00 1 1h
Number of Iterations....: 15
(scaled) (unscaled)
Objective...............: 4.7952701409621918e+04 4.7952701409621918e+04
Dual infeasibility......: 8.1637622585163618e-19 8.1637622585163618e-19
Constraint violation....: 2.3083090823202061e-21 2.3083090823202061e-21
Complementarity.........: 1.1469774205081724e-18 1.1469774205081724e-18
Overall NLP error.......: 1.1469774205081724e-18 1.1469774205081724e-18
Number of objective function evaluations = 16
Number of objective gradient evaluations = 16
Number of constraint evaluations = 16
Number of constraint Jacobian evaluations = 16
Number of Lagrangian Hessian evaluations = 15
Number of KKT factorizations = 16
Number of KKT backsolves = 16
Total wall secs in initialization = 4.840 s
Total wall secs in linear solver = 0.576 s
Total wall secs in NLP function evaluations = 0.034 s
Total wall secs in solver (w/o init./fun./lin. alg.) = 19.161 s
Total wall secs = 24.611 s
EXIT: Optimal Solution Found (tol = 1.0e-17).Note that the final tolerance is much lower than before. We get the solution in quadruple precision
results_128.solution84-element Vector{Quadmath.Float128}:
-6.10422346779984075324746317637300190e+00
-7.32016751646112951214547122344744191e+00
-8.70225021519367902455800396043161132e+00
-6.54778300598681716679558184769648358e+00
-5.15677166066576857882792575858440579e+00
-1.84695048479262641750440403327346213e+00
-3.32537558953529983731321975183318410e+00
-4.92515250359083516678377485953970913e-01
1.35832327690916348435514663903931195e+00
3.82340248568206791959202066079244177e+00
⋮
-2.68405535104923868206466270293099121e+00
-3.14226949616395173158877302605268355e+00
-4.00022305007620456953394893097651818e+00
-4.97642231990682172164846696875993814e+00
-6.14153113112896516196594293898607652e+00
-7.10621611493401363964299609878440189e+00
-5.74357155980421850420980936002355446e+00
-4.61682282079019122893368919072475892e+00
3.09035786082032038293382403877706561e-01as well as the final objective:
results_128.objective4.79527014096219181130758389630855969e+04