[API] Expose the options for iterative refinements and quasi-Newton

src/IPM/IPM.jl

@@ -17,7 +17,7 @@ mutable struct MadNLPSolver{
     IC <: AbstractInertiaCorrector,
     KKTVec <: AbstractKKTVector{T, VT}
     } <: AbstractMadNLPSolver{T}
-    
+
     nlp::Model
     cb::CB
     kkt::KKTSystem
@@ -103,65 +103,66 @@ mutable struct MadNLPSolver{
 end
 
 function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT}
-
-    opt, opt_linear_solver, logger = load_options(nlp; kwargs...)
-    @assert is_supported(opt.linear_solver, T)
+    options = load_options(nlp; kwargs...)
+
+    ipm_opt = options.interior_point
+    @assert is_supported(ipm_opt.linear_solver, T)
 
     cnt = MadNLPCounters(start_time=time())
-    cb = create_callback(opt.callback, nlp, opt)
-    
+    cb = create_callback(ipm_opt.callback, nlp, ipm_opt)
+
     # generic options
-    opt.disable_garbage_collector &&
+    ipm_opt.disable_garbage_collector &&
         (GC.enable(false); @warn(logger,"Julia garbage collector is temporarily disabled"))
-    set_blas_num_threads(opt.blas_num_threads; permanent=true)
-    @trace(logger,"Initializing variables.")
-    
+    set_blas_num_threads(ipm_opt.blas_num_threads; permanent=true)
+    @trace(options.logger,"Initializing variables.")
+
     ind_cons = get_index_constraints(
         get_lvar(nlp), get_uvar(nlp),
         get_lcon(nlp), get_ucon(nlp),
-        opt.fixed_variable_treatment,
-        opt.equality_treatment
+        ipm_opt.fixed_variable_treatment,
+        ipm_opt.equality_treatment
     )
 
     ind_lb = ind_cons.ind_lb
     ind_ub = ind_cons.ind_ub
-    
+
     ns = length(ind_cons.ind_ineq)
     nx = get_nvar(nlp)
     n = nx+ns
     m = get_ncon(nlp)
     nlb = length(ind_lb)
     nub = length(ind_ub)
 
-    @trace(logger,"Initializing KKT system.")
+    @trace(options.logger,"Initializing KKT system.")
     kkt = create_kkt_system(
-        opt.kkt_system,
+        ipm_opt.kkt_system,
         cb,
-        opt,
-        opt_linear_solver,
+        ipm_opt,
+        options.linear_solver,
         cnt,
         ind_cons
     )
 
-    @trace(logger,"Initializing iterative solver.")
-    iterator = opt.iterator(kkt; cnt = cnt, logger = logger)
+    @trace(options.logger,"Initializing iterative solver.")
+    iterator = ipm_opt.iterator(kkt; cnt = cnt, logger = options.logger, opt = options.iterative_refinement)
 
     x = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
     xl = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
-    xu = PrimalVector(VT, nx, ns, ind_lb, ind_ub)  
+    xu = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
     zl = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
     zu = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
     f = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
     x_trial = PrimalVector(VT, nx, ns, ind_lb, ind_ub)
-    
+
     d = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
     p = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
     _w1 = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
     _w2 = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
     _w3 = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
     _w4 = UnreducedKKTVector(VT, n, m, nlb, nub, ind_lb, ind_ub)
 
-    jacl = VT(undef,n) 
+    jacl = VT(undef,n)
     c_trial = VT(undef, m)
     y = VT(undef, m)
     c = VT(undef, m)
@@ -179,39 +180,39 @@ function MadNLPSolver(nlp::AbstractNLPModel{T,VT}; kwargs...) where {T, VT}
     dx_lr = view(d.xp, ind_cons.ind_lb) # TODO
     dx_ur = view(d.xp, ind_cons.ind_ub) # TODO
 
-    inertia_correction_method = if opt.inertia_correction_method == InertiaAuto
+    inertia_correction_method = if ipm_opt.inertia_correction_method == InertiaAuto
         is_inertia(kkt.linear_solver)::Bool ? InertiaBased : InertiaFree
     else
-        opt.inertia_correction_method
+        ipm_opt.inertia_correction_method
     end
 
     inertia_corrector = build_inertia_corrector(
         inertia_correction_method,
         VT,
         n, m, nlb, nub, ind_lb, ind_ub
     )
-    
+
     cnt.init_time = time() - cnt.start_time
 
     return MadNLPSolver(
         nlp, cb, kkt,
-        opt, cnt, logger, 
+        ipm_opt, cnt, options.logger,
         n, m, nlb, nub,
         x, y, zl, zu, xl, xu,
-        zero(T), f, c, 
-        jacl, 
-        d, p, 
-        _w1, _w2, _w3, _w4, 
-        x_trial, c_trial, zero(T), c_slk, rhs, 
-        ind_cons.ind_ineq, ind_cons.ind_fixed, ind_cons.ind_llb, ind_cons.ind_uub, 
-        x_lr, x_ur, xl_r, xu_r, zl_r, zu_r, dx_lr, dx_ur, x_trial_lr, x_trial_ur, 
-        iterator, 
+        zero(T), f, c,
+        jacl,
+        d, p,
+        _w1, _w2, _w3, _w4,
+        x_trial, c_trial, zero(T), c_slk, rhs,
+        ind_cons.ind_ineq, ind_cons.ind_fixed, ind_cons.ind_llb, ind_cons.ind_uub,
+        x_lr, x_ur, xl_r, xu_r, zl_r, zu_r, dx_lr, dx_ur, x_trial_lr, x_trial_ur,
+        iterator,
         zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T), zero(T),
         " ",
         zero(T), zero(T), zero(T),
         Tuple{T, T}[],
         inertia_corrector, nothing,
-        INITIAL, Dict(), 
+        INITIAL, Dict(),
     )
 
 end

src/KKT/dense.jl

@@ -76,7 +76,10 @@ function create_kkt_system(
     fill!(du_diag, zero(T))
     fill!(diag_hess, zero(T))
 
-    quasi_newton = create_quasi_newton(opt.hessian_approximation, cb, n)
+    quasi_newton = create_quasi_newton(
+        opt.hessian_approximation, cb, n;
+        options=opt.quasi_newton_options,
+    )
     cnt.linear_solver_time += @elapsed linear_solver = opt.linear_solver(aug_com; opt = opt_linear_solver)
 
     return DenseKKTSystem(
@@ -179,7 +182,10 @@ function create_kkt_system(
     ind_eq_shifted = ind_cons.ind_eq .+ n .+ ns
     ind_ineq_shifted = ind_cons.ind_ineq .+ n .+ ns
 
-    quasi_newton = create_quasi_newton(opt.hessian_approximation, cb, n)
+    quasi_newton = create_quasi_newton(
+        opt.hessian_approximation, cb, n;
+        options=opt.quasi_newton_options,
+    )
     cnt.linear_solver_time += @elapsed linear_solver = opt.linear_solver(aug_com; opt = opt_linear_solver)
 
     return DenseCondensedKKTSystem(

src/KKT/sparse.jl

@@ -206,7 +206,10 @@ function create_kkt_system(
     jac_sparsity_J = create_array(cb, Int32, cb.nnzj)
     _jac_sparsity_wrapper!(cb,jac_sparsity_I, jac_sparsity_J)
 
-    quasi_newton = create_quasi_newton(opt.hessian_approximation, cb, n)
+    quasi_newton = create_quasi_newton(
+        opt.hessian_approximation, cb, n;
+        options=opt.quasi_newton_options,
+    )
     hess_sparsity_I, hess_sparsity_J = build_hessian_structure(cb, opt.hessian_approximation)
 
     nlb = length(ind_cons.ind_lb)
@@ -316,7 +319,10 @@ function create_kkt_system(
     m = cb.ncon
 
     # Quasi-newton
-    quasi_newton = create_quasi_newton(opt.hessian_approximation, cb, n)
+    quasi_newton = create_quasi_newton(
+        opt.hessian_approximation, cb, n;
+        options=opt.quasi_newton_options,
+    )
 
     # Evaluate sparsity pattern
     jac_sparsity_I = create_array(cb, Int32, cb.nnzj)
@@ -474,7 +480,10 @@ function create_kkt_system(
     jac_sparsity_J = create_array(cb, Int32, cb.nnzj)
     _jac_sparsity_wrapper!(cb,jac_sparsity_I, jac_sparsity_J)
 
-    quasi_newton = create_quasi_newton(opt.hessian_approximation, cb, n)
+    quasi_newton = create_quasi_newton(
+        opt.hessian_approximation, cb, n;
+        options=opt.quasi_newton_options,
+    )
     hess_sparsity_I, hess_sparsity_J = build_hessian_structure(cb, opt.hessian_approximation)
 
     force_lower_triangular!(hess_sparsity_I,hess_sparsity_J)

src/LinearSolvers/backsolve.jl

@@ -1,7 +1,7 @@
-@kwdef struct RichardsonOptions
-    max_iter::Int = 10
-    tol::Float64 = 1e-10
-    acceptable_tol::Float64 = 1e-5
+@kwdef mutable struct RichardsonOptions <: AbstractOptions
+    richardson_max_iter::Int = 10
+    richardson_tol::Float64 = 1e-10
+    richardson_acceptable_tol::Float64 = 1e-5
 end
 
 struct RichardsonIterator{T, KKT <: AbstractKKTSystem{T}} <: AbstractIterator{T}
@@ -22,6 +22,8 @@ function RichardsonIterator(
     )
 end
 
+default_options(::Type{RichardsonIterator}) = RichardsonOptions()
+
 function solve_refine!(
     x::VT,
     iterator::R,
@@ -66,7 +68,7 @@ function solve_refine!(
         @debug(iterator.logger, @sprintf("%4i %6.2e", iter, residual_ratio))
         iter += 1
 
-        if (iter >= iterator.opt.max_iter) || (residual_ratio < iterator.opt.tol)
+        if (iter >= iterator.opt.richardson_max_iter) || (residual_ratio < iterator.opt.richardson_tol)
             break
         end
     end
@@ -79,6 +81,6 @@ function solve_refine!(
         ),
     )
 
-    return residual_ratio < iterator.opt.acceptable_tol
+    return residual_ratio < iterator.opt.richardson_acceptable_tol
 end
 

src/options.jl

@@ -47,6 +47,7 @@ end
     hessian_constant::Bool = false
     kkt_system::Type = SparseKKTSystem
     hessian_approximation::Type = ExactHessian
+    quasi_newton_options::QuasiNewtonOptions = QuasiNewtonOptions()
     callback::Type = SparseCallback
 
     # initialization options
@@ -163,7 +164,10 @@ function load_options(nlp; options...)
     check_option_sanity(opt_ipm)
     # Initiate linear-solver options
     opt_linear_solver = default_options(opt_ipm.linear_solver)
-    remaining_options = set_options!(opt_linear_solver, linear_solver_options)
+    iterator_options = set_options!(opt_linear_solver, linear_solver_options)
+    # Initiate iterator options
+    opt_iterator = default_options(opt_ipm.iterator)
+    remaining_options = set_options!(opt_iterator, iterator_options)
 
     # Initiate logger
     logger = MadNLPLogger(
@@ -177,6 +181,11 @@ function load_options(nlp; options...)
     if !isempty(remaining_options)
         print_ignored_options(logger, remaining_options)
     end
-    return opt_ipm, opt_linear_solver, logger
+    return (
+        interior_point=opt_ipm,
+        linear_solver=opt_linear_solver,
+        iterative_refinement=opt_iterator,
+        logger=logger,
+    )
 end
 

src/quasi_newton.jl

@@ -35,6 +35,11 @@ curvature(::Val{SCALAR2}, sk, yk) = dot(yk, yk) / dot(sk, yk)
 curvature(::Val{SCALAR3}, sk, yk) = 0.5 * (curvature(Val(SCALAR1), sk, yk) + curvature(Val(SCALAR2), sk, yk))
 curvature(::Val{SCALAR4}, sk, yk) = sqrt(curvature(Val(SCALAR1), sk, yk) * curvature(Val(SCALAR2), sk, yk))
 
+@kwdef mutable struct QuasiNewtonOptions <: AbstractOptions
+    init_strategy::BFGSInitStrategy = SCALAR1
+    max_history::Int = 6
+end
+
 
 """
     BFGS{T, VT} <: AbstractQuasiNewton{T, VT}
@@ -61,10 +66,10 @@ function create_quasi_newton(
     ::Type{BFGS},
     cb::AbstractCallback{T,VT},
     n;
-    init_strategy = SCALAR1
+    options=QuasiNewtonOptions(),
     ) where {T,VT}
     BFGS(
-        init_strategy,
+        options.init_strategy,
         VT(undef, n),
         VT(undef, n),
         VT(undef, n),
@@ -100,10 +105,10 @@ function create_quasi_newton(
     ::Type{DampedBFGS},
     cb::AbstractCallback{T,VT},
     n;
-    init_strategy = SCALAR1
+    options=QuasiNewtonOptions(),
     ) where {T,VT}
     return DampedBFGS(
-        init_strategy,
+        options.init_strategy,
         VT(undef, n),
         VT(undef, n),
         VT(undef, n),
@@ -178,17 +183,16 @@ function create_quasi_newton(
     ::Type{CompactLBFGS},
     cb::AbstractCallback{T,VT},
     n;
-    max_mem=6,
-    init_strategy = SCALAR1
+    options=QuasiNewtonOptions(),
     ) where {T, VT}
     return CompactLBFGS(
-        init_strategy,
+        options.init_strategy,
         fill!(create_array(cb, n), zero(T)),
         fill!(create_array(cb, n), zero(T)),
         fill!(create_array(cb, n), zero(T)),
         fill!(create_array(cb, n), zero(T)),
         fill!(create_array(cb, n), zero(T)),
-        max_mem,
+        options.max_history,
         0,
         fill!(create_array(cb, n, 0), zero(T)),
         fill!(create_array(cb, n, 0), zero(T)),
@@ -319,4 +323,4 @@ end
 
 
 struct ExactHessian{T, VT} <: AbstractHessian{T, VT} end
-create_quasi_newton(::Type{ExactHessian}, cb::AbstractCallback{T,VT}, n) where {T,VT} = ExactHessian{T, VT}()
+create_quasi_newton(::Type{ExactHessian}, cb::AbstractCallback{T,VT}, n; options...) where {T,VT} = ExactHessian{T, VT}()