diff --git a/Project.toml b/Project.toml index d5005d1a..d195f3c3 100644 --- a/Project.toml +++ b/Project.toml @@ -1,11 +1,6 @@ name = "RegularizedOptimization" uuid = "20620ad1-4fe4-4467-ae46-fb087718fe7b" -author = ["Robert Baraldi ", - "Youssef Diouane ", - "Maxence Gollier ", - "Mohamed Laghdaf Habiboullah ", - "Geoffroy Leconte ", - "Dominique Orban "] +author = ["Robert Baraldi ", "Youssef Diouane ", "Maxence Gollier ", "Mohamed Laghdaf Habiboullah ", "Geoffroy Leconte ", "Dominique Orban "] version = "0.1.0" [deps] @@ -18,6 +13,7 @@ NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProxTV = "925ea013-038b-5ab6-a1ab-e0849925e528" ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537" RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278" ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263" @@ -32,6 +28,7 @@ NLPModels = "0.19, 0.20, 0.21" NLPModelsModifiers = "0.7" OptimizationProblems = "0.9.2" Percival = "0.7.2" +ProxTV = "1.2.3" ProximalOperators = "0.15" RegularizedProblems = "0.1.1" ShiftedProximalOperators = "0.2" diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl index 85781388..d91063d1 100644 --- a/src/RegularizedOptimization.jl +++ b/src/RegularizedOptimization.jl @@ -13,7 +13,8 @@ using LinearOperators, NLPModelsModifiers, RegularizedProblems, ShiftedProximalOperators, - SolverCore + SolverCore, + ProxTV using Percival: AugLagModel, update_y!, update_μ! import SolverCore.reset! @@ -46,7 +47,9 @@ include("LM_alg.jl") include("LMTR_alg.jl") include("R2DH.jl") include("R2NModel.jl") +include("iR2_alg.jl") include("R2N.jl") include("AL_alg.jl") +include("iR2N.jl") end # module RegularizedOptimization diff --git a/src/iR2N.jl b/src/iR2N.jl new file mode 100644 index 00000000..30a9ecc9 --- /dev/null +++ b/src/iR2N.jl @@ -0,0 +1,596 @@ +export iR2N, iR2NSolver, solve! + +import SolverCore.solve! + +mutable struct iR2NSolver{ + T <: Real, + G <: Union{ShiftedProximableFunction, InexactShiftedProximableFunction}, + V <: AbstractVector{T}, + ST <: AbstractOptimizationSolver, + PB <: AbstractRegularizedNLPModel, +} <: AbstractOptimizationSolver + xk::V + ∇fk::V + ∇fk⁻::V + mν∇fk::V + ψ::G + xkn::V + s::V + s1::V + has_bnds::Bool + l_bound::V + u_bound::V + l_bound_m_x::V + u_bound_m_x::V + m_fh_hist::V + subsolver::ST + subpb::PB + substats::GenericExecutionStats{T, V, V, T} +end + +function iR2NSolver( + reg_nlp::AbstractRegularizedNLPModel{T, V}; + subsolver = iR2Solver, + m_monotone::Int = 1, +) where {T, V} + x0 = reg_nlp.model.meta.x0 + l_bound = reg_nlp.model.meta.lvar + u_bound = reg_nlp.model.meta.uvar + + xk = similar(x0) + ∇fk = similar(x0) + ∇fk⁻ = similar(x0) + mν∇fk = similar(x0) + xkn = similar(x0) + s = similar(x0) + s1 = similar(x0) + v = similar(x0) + has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf)) + if has_bnds + l_bound_m_x = similar(xk) + u_bound_m_x = similar(xk) + @. l_bound_m_x = l_bound - x0 + @. u_bound_m_x = u_bound - x0 + else + l_bound_m_x = similar(xk, 0) + u_bound_m_x = similar(xk, 0) + end + m_fh_hist = fill(T(-Inf), m_monotone - 1) + + ψ = + has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : + shifted(reg_nlp.h, xk) + + if ψ isa InexactShiftedProximableFunction + + # Create a completely new context for the subsolver + function create_new_context(h, xk, l_bound_m_x, u_bound_m_x, selected) + new_h = deepcopy(h) + new_context = deepcopy(h.context) + new_h.context = new_context + return has_bnds ? shifted(new_h, xk, l_bound_m_x, u_bound_m_x, selected) : shifted(new_h, xk) + end + + ψ_sub = create_new_context(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) + else + ψ_sub = ψ + end + + Bk = hess_op(reg_nlp.model, x0) + sub_nlp = R2NModel(Bk, ∇fk, T(1), x0) + subpb = RegularizedNLPModel(sub_nlp, ψ_sub) + substats = RegularizedExecutionStats(subpb) + subsolver = subsolver(subpb) + + return iR2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subpb)}( + xk, + ∇fk, + ∇fk⁻, + mν∇fk, + ψ, + xkn, + s, + s1, + has_bnds, + l_bound, + u_bound, + l_bound_m_x, + u_bound_m_x, + m_fh_hist, + subsolver, + subpb, + substats, + ) +end + +""" + iR2N(reg_nlp; kwargs…) + +A second-order quadratic regularization method for the problem + + min f(x) + h(x) + +where f: ℝⁿ → ℝ is C¹, and h: ℝⁿ → ℝ is +lower semi-continuous, proper and prox-bounded. + +About each iterate xₖ, a step sₖ is computed as an approximate solution of + + min φ(s; xₖ) + ½ σₖ ‖s‖² + ψ(s; xₖ) + +where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs + ½ sᵀBₖs is a quadratic approximation of f about xₖ, +ψ(s; xₖ) is either h(xₖ + s) or an approximation of h(xₖ + s), ‖⋅‖ is the ℓ₂ norm and σₖ > 0 is the regularization parameter. + +For advanced usage, first define a solver "iR2NSolver" to preallocate the memory used in the algorithm, and then call `solve!`: + + solver = iR2NSolver(reg_nlp; subsolver = iR2Solver, m_monotone = 1) + solve!(solver, reg_nlp) + + stats = RegularizedExecutionStats(reg_nlp) + solve!(solver, reg_nlp, stats) + +# Arguments +* `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the problem to solve, see `RegularizedProblems.jl`, `NLPModels.jl`. + +# Keyword arguments +- `x::V = nlp.meta.x0`: the initial guess; +- `atol::T = √eps(T)`: absolute tolerance; +- `rtol::T = √eps(T)`: relative tolerance; +- `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance +- `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); +- `max_time::Float64 = 30.0`: maximum time limit in seconds; +- `max_iter::Int = 10000`: maximum number of iterations; +- `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration; +- `σmin::T = eps(T)`: minimum value of the regularization parameter; +- `η1::T = √√eps(T)`: successful iteration threshold; +- `η2::T = T(0.9)`: very successful iteration threshold; +- `ν::T = eps(T)^(1 / 5)`: inverse of the initial regularization parameter: ν = 1/σ; +- `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful. +- `θ::T = 1/(1 + eps(T)^(1 / 5))`: is the model decrease fraction with respect to the decrease of the Cauchy model. +- `m_monotone::Int = 1`: monotonicity parameter. By default, iR2N is monotone but the non-monotone variant will be used if `m_monotone > 1` + +The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure. + +# Output +The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. + +# Callback +$(callback_docstring) +""" +function iR2N( + nlp::AbstractNLPModel{T, V}, + h, + options::ROSolverOptions{T}; + kwargs..., +) where {T <: Real, V} + kwargs_dict = Dict(kwargs...) + selected = pop!(kwargs_dict, :selected, 1:(nlp.meta.nvar)) + x0 = pop!(kwargs_dict, :x0, nlp.meta.x0) + reg_nlp = RegularizedNLPModel(nlp, h, selected) + return iR2N( + reg_nlp, + x = x0, + atol = options.ϵa, + rtol = options.ϵr, + neg_tol = options.neg_tol, + verbose = options.verbose, + max_iter = options.maxIter, + max_time = options.maxTime, + σmin = options.σmin, + η1 = options.η1, + η2 = options.η2, + ν = options.ν, + γ = options.γ; + kwargs_dict..., + ) +end + +function iR2N(reg_nlp::AbstractRegularizedNLPModel; kwargs...) + if !(shifted(reg_nlp.h, reg_nlp.model.meta.x0) isa InexactShiftedProximableFunction) + @warn "h has exact prox, switching to R2N" + return R2N(reg_nlp; kwargs...) + end + kwargs_dict = Dict(kwargs...) + m_monotone = pop!(kwargs_dict, :m_monotone, 1) + subsolver = pop!(kwargs_dict, :subsolver, iR2Solver) + solver = iR2NSolver(reg_nlp, subsolver = subsolver, m_monotone = m_monotone) + stats = GenericExecutionStats(reg_nlp.model) + solve!(solver, reg_nlp, stats; kwargs_dict...) + return stats +end + +function SolverCore.solve!( + solver::iR2NSolver{T, G, V}, + reg_nlp::AbstractRegularizedNLPModel{T, V}, + stats::GenericExecutionStats{T, V}; + callback = (args...) -> nothing, + x::V = reg_nlp.model.meta.x0, + atol::T = √eps(T), + rtol::T = √eps(T), + neg_tol::T = eps(T)^(1 / 4), + verbose::Int = 0, + max_iter::Int = 10000, + max_time::Float64 = 30.0, + max_eval::Int = -1, + σmin::T = eps(T), + η1::T = √√eps(T), + η2::T = T(0.9), + ν::T = eps(T)^(1 / 5), + γ::T = T(3), + β::T = 1 / eps(T), + θ::T = 1 / (1 + eps(T)^(1 / 5)), + kwargs..., +) where {T, V, G} + reset!(stats) + + # Retrieve workspace + selected = reg_nlp.selected + h = reg_nlp.h + nlp = reg_nlp.model + + xk = solver.xk .= x + + # Make sure ψ has the correct shift + shift!(solver.ψ, xk) + + ∇fk = solver.∇fk + ∇fk⁻ = solver.∇fk⁻ + mν∇fk = solver.mν∇fk + ψ = solver.ψ + if ψ isa ShiftedProximableFunction + κs = 1.0 + else + κs = ψ.h.context.κs + end + xkn = solver.xkn + s = solver.s + s1 = solver.s1 + m_fh_hist = solver.m_fh_hist .= T(-Inf) + has_bnds = solver.has_bnds + + if has_bnds + l_bound_m_x = solver.l_bound_m_x + u_bound_m_x = solver.u_bound_m_x + l_bound = solver.l_bound + u_bound = solver.u_bound + end + m_monotone = length(m_fh_hist) + 1 + + # initialize parameters + improper = false + hk = @views h(xk[selected]) + if hk == Inf + verbose > 0 && @info "iR2N: finding initial guess where nonsmooth term is finite" + prox!(xk, h, xk, T(1)) + hk = @views h(xk[selected]) + hk < Inf || error("prox computation must be erroneous") + verbose > 0 && @debug "iR2N: found point where h has value" hk + end + improper = (hk == -Inf) + improper == true && @warn "iR2N: Improper term detected" + improper == true && return stats + + if verbose > 0 + @info log_header( + [:outer, :inner, :fx, :hx, :norm_s1_νInv, :ρ, :σ, :normx, :norms, :normB, :arrow], + [Int, Int, T, T, T, T, T, T, T, T, Char], + hdr_override = Dict{Symbol, String}( + :fx => "f(x)", + :hx => "h(x)", + :norm_s1_νInv => "‖s‖/ν", + :normx => "‖x‖", + :norms => "‖s‖", + :normB => "‖B‖", + :arrow => "iR2N", + ), + colsep = 1, + ) + end + + local ξ1::T + local ρk::T = zero(T) + + σk = max(1 / ν, σmin) + + fk = obj(nlp, xk) + grad!(nlp, xk, solver.∇fk) + ∇fk⁻ .= solver.∇fk + + quasiNewtTest = isa(nlp, QuasiNewtonModel) + λmax::T = T(1) + solver.subpb.model.B = hess_op(nlp, xk) + + λmax, found_λ = opnorm(solver.subpb.model.B) + found_λ || error("operator norm computation failed") + + ν₁ = θ / (λmax + σk) + ν_sub = ν₁ + + sqrt_ξ1_νInv = one(T) + + @. mν∇fk = -ν₁ * solver.∇fk + + set_iter!(stats, 0) + start_time = time() + set_time!(stats, 0.0) + set_objective!(stats, fk + hk) + set_solver_specific!(stats, :smooth_obj, fk) + set_solver_specific!(stats, :nonsmooth_obj, hk) + m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk) + + φ1 = let ∇fk = solver.∇fk + d -> dot(∇fk, d) + end + + mk1 = let ψ = ψ + d -> φ1(d) + ψ(d)::T + end + + mk = let ψ = ψ, solver = solver + d -> obj(solver.subpb.model, d) + ψ(d)::T + end + + update_prox_context!(solver, stats, ψ) + prox!(s1, ψ, mν∇fk, ν₁) + mks = mk1(s1) + + # Record outer-level prox call counters (exclude feasibility prox above) + if ψ isa InexactShiftedProximableFunction + ctx = ψ.h.context + ctx.prox_calls_outer += 1 + ctx.prox_iters_outer += ctx.last_prox_iters + end + + ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps() + sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) + norm_s1_νInv = norm(s1) / ν₁ + # solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol * (√κs)) + solved = norm_s1_νInv ≤ atol * κs # new stopping test according to theory + (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && + @warn "iR2N: prox-gradient step should produce a decrease but ξ1 = $(ξ1) and sqrt_ξ1_νInv = $sqrt_ξ1_νInv > $neg_tol" + atol += rtol * norm_s1_νInv # make stopping test absolute and relative + + set_status!( + stats, + get_status( + reg_nlp, + elapsed_time = stats.elapsed_time, + iter = stats.iter, + optimal = solved, + improper = improper, + max_eval = max_eval, + max_time = max_time, + max_iter = max_iter, + ), + ) + + callback(nlp, solver, stats) + + done = stats.status != :unknown + + while !done + sub_atol = stats.iter == 0 ? 1.0e-3 : min(norm_s1_νInv^(1.5), norm_s1_νInv * 1e-3) + + solver.subpb.model.σ = σk + isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1 / ν₁) + ν_sub = isa(solver.subsolver, R2DHSolver) ? 1 / σk : ν₁ + # Snapshot inner context counters to aggregate deltas only + inner_ctx = solver.subpb.h.h.context + prev_calls_inner = inner_ctx.prox_calls_inner + prev_iters_inner = inner_ctx.prox_iters_inner + prev_proxstats3_inner = inner_ctx.prox_stats[3] + + solve!( + solver.subsolver, + solver.subpb, + solver.substats; + x = s1, + atol = sub_atol, + ν = ν_sub, + neg_tol = neg_tol, + kwargs..., + ) + + s .= solver.substats.solution + + if ψ isa InexactShiftedProximableFunction + # Aggregate inner-level prox counters by delta since last call + calls_delta = inner_ctx.prox_calls_inner - prev_calls_inner + iters_delta = inner_ctx.prox_iters_inner - prev_iters_inner + iters_prox_delta = inner_ctx.prox_stats[3] - prev_proxstats3_inner + let outer_ctx = ψ.h.context + outer_ctx.prox_calls_inner += calls_delta + outer_ctx.prox_iters_inner += iters_delta + outer_ctx.prox_stats[3] += iters_prox_delta + # Keep subsolver iteration total for legacy metric + outer_ctx.prox_stats[2] += solver.substats.iter + end + end + + if norm(s) > β * norm(s1) + s .= s1 + end + + xkn .= xk .+ s + fkn = obj(nlp, xkn) + hkn = @views h(xkn[selected]) + mks = mk(s) + + fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk + Δobj = fhmax - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps() + Δmod = fhmax - (fk + mks) + max(1, abs(fhmax)) * 10 * eps() + ξ = hk - mks + max(1, abs(hk)) * 10 * eps() + sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν₁) : sqrt(-ξ / ν₁) + + if (ξ < 0 && sqrt_ξ_νInv > neg_tol) || isnan(ξ) + error("iR2N: failed to compute a step: ξ = $ξ and sqrt_ξ_νInv = $sqrt_ξ_νInv") + end + + aux_assert = (1 - θ) * ξ1 + if ξ < aux_assert && sqrt_ξ_νInv > neg_tol + @warn("iR2N: ξ should be ≥ (1 - θ) * ξ1 but ξ = $ξ and (1 - θ) * ξ1 = $aux_assert.",) + end + ρk = Δobj / Δmod + + verbose > 0 && + stats.iter % verbose == 0 && + @info log_row( + Any[ + stats.iter, + solver.substats.iter, + fk, + hk, + norm_s1_νInv, + ρk, + σk, + norm(xk), + norm(s), + λmax, + (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "="), + ], + colsep = 1, + ) + + if η1 ≤ ρk < Inf + xk .= xkn + if has_bnds + @. l_bound_m_x = l_bound - xk + @. u_bound_m_x = u_bound - xk + set_bounds!(ψ, l_bound_m_x, u_bound_m_x) + end + #update functions + fk = fkn + hk = hkn + + shift!(ψ, xk) + grad!(nlp, xk, solver.∇fk) + + if quasiNewtTest + @. ∇fk⁻ = solver.∇fk - ∇fk⁻ + push!(nlp, s, ∇fk⁻) + end + solver.subpb.model.B = hess_op(nlp, xk) + + λmax, found_λ = opnorm(solver.subpb.model.B) + found_λ || error("operator norm computation failed") + + ∇fk⁻ .= solver.∇fk + end + + if η2 ≤ ρk < Inf + σk = max(σk / γ, σmin) + end + + if ρk < η1 || ρk == Inf + σk = σk * γ + end + + ν₁ = θ / (λmax + σk) + m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk) + + set_objective!(stats, fk + hk) + set_solver_specific!(stats, :smooth_obj, fk) + set_solver_specific!(stats, :nonsmooth_obj, hk) + set_iter!(stats, stats.iter + 1) + set_time!(stats, time() - start_time) + + @. mν∇fk = -ν₁ * solver.∇fk + + update_prox_context!(solver, stats, ψ) + prox!(s1, ψ, mν∇fk, ν₁) + mks = mk1(s1) + + # Record outer-level prox call counters per iteration + if ψ isa InexactShiftedProximableFunction + ctx = ψ.h.context + ctx.prox_calls_outer += 1 + ctx.prox_iters_outer += ctx.last_prox_iters + end + + ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps() + sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) + norm_s1_νInv = norm(s1) / ν₁ + # solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol * (√κs)) + solved = norm_s1_νInv ≤ atol * κs # new stopping test according to theory + + (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && + @warn "iR2N: prox-gradient step should produce a decrease but ξ1 = $(ξ1) with sqrt_ξ1_νInv = $sqrt_ξ1_νInv > $neg_tol" + set_status!( + stats, + get_status( + reg_nlp, + elapsed_time = stats.elapsed_time, + iter = stats.iter, + optimal = solved, + improper = improper, + max_eval = max_eval, + max_time = max_time, + max_iter = max_iter, + ), + ) + + callback(nlp, solver, stats) + + done = stats.status != :unknown + end + + if verbose > 0 && stats.status == :first_order + @info log_row( + Any[ + stats.iter, + 0, + fk, + hk, + norm_s1_νInv, + ρk, + σk, + norm(xk), + norm(s), + λmax, + (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "="), + ], + colsep = 1, + ) + @info "iR2N: terminating with ‖s‖/ν = $(norm_s1_νInv)" + end + + if ψ isa InexactShiftedProximableFunction + set_solver_specific!(stats, :total_iters_subsolver, solver.ψ.h.context.prox_stats[2]) + set_solver_specific!( + stats, + :mean_iters_subsolver, + solver.ψ.h.context.prox_stats[2] / stats.iter, + ) + set_solver_specific!(stats, :total_iters_prox, solver.ψ.h.context.prox_stats[3]) + set_solver_specific!(stats, :mean_iters_prox, solver.ψ.h.context.prox_stats[3] / stats.iter) #TODO: work on this line to specify iR2 prox calls and iR2N prox calls + + # Detailed outer/inner/all prox stats (per-call means) + let ctx = solver.ψ.h.context + total_calls_outer = ctx.prox_calls_outer + total_iters_outer = ctx.prox_iters_outer + mean_outer = total_calls_outer > 0 ? total_iters_outer / total_calls_outer : zero(eltype(xk)) + + total_calls_inner = ctx.prox_calls_inner + total_iters_inner = ctx.prox_iters_inner + mean_inner = total_calls_inner > 0 ? total_iters_inner / total_calls_inner : zero(eltype(xk)) + + total_calls_all = total_calls_outer + total_calls_inner + total_iters_all = total_iters_outer + total_iters_inner + mean_all = total_calls_all > 0 ? total_iters_all / total_calls_all : zero(eltype(xk)) + + set_solver_specific!(stats, :prox_calls_outer, total_calls_outer) + set_solver_specific!(stats, :total_iters_prox_outer, total_iters_outer) + set_solver_specific!(stats, :mean_iters_prox_outer, mean_outer) + + set_solver_specific!(stats, :prox_calls_inner, total_calls_inner) + set_solver_specific!(stats, :total_iters_prox_inner, total_iters_inner) + set_solver_specific!(stats, :mean_iters_prox_inner, mean_inner) + + set_solver_specific!(stats, :prox_calls_all, total_calls_all) + set_solver_specific!(stats, :total_iters_prox_all, total_iters_all) + set_solver_specific!(stats, :mean_iters_prox_all, mean_all) + end + end + set_solution!(stats, xk) + set_residuals!(stats, zero(eltype(xk)), norm_s1_νInv) + return stats +end \ No newline at end of file diff --git a/src/iR2_alg.jl b/src/iR2_alg.jl new file mode 100644 index 00000000..7d098d3b --- /dev/null +++ b/src/iR2_alg.jl @@ -0,0 +1,617 @@ +export iR2, iR2Solver, solve! + +import SolverCore.solve! + +mutable struct iR2Solver{ + R <: Real, + G <: Union{ShiftedProximableFunction, InexactShiftedProximableFunction, Nothing}, + S <: AbstractVector{R}, +} <: AbstractOptimizationSolver + xk::S + ∇fk::S + mν∇fk::S + ψ::G + xkn::S + s::S + has_bnds::Bool + l_bound::S + u_bound::S + l_bound_m_x::S + u_bound_m_x::S + Fobj_hist::Vector{R} + Hobj_hist::Vector{R} + Complex_hist::Vector{Int} +end + +# !!!!!! Not used anywhere !!!!! +# function iR2Solver( +# x0::S, +# options::ROSolverOptions, +# l_bound::S, +# u_bound::S; +# ψ = nothing, +# ) where {R <: Real, S <: AbstractVector{R}} +# maxIter = options.maxIter +# xk = similar(x0) +# ∇fk = similar(x0) +# mν∇fk = similar(x0) +# xkn = similar(x0) +# s = zero(x0) +# has_bnds = any(l_bound .!= R(-Inf)) || any(u_bound .!= R(Inf)) +# if has_bnds +# l_bound_m_x = similar(xk) +# u_bound_m_x = similar(xk) +# else +# l_bound_m_x = similar(xk, 0) +# u_bound_m_x = similar(xk, 0) +# end +# Fobj_hist = zeros(R, maxIter + 2) +# Hobj_hist = zeros(R, maxIter + 2) +# Complex_hist = zeros(Int, maxIter + 2) +# dualGap = options.dualGap +# κs = options.κs +# mk1 = options.mk1 +# shift = similar(x0) +# s_k_unshifted = similar(x0) +# callback_pointer = options.callback_pointer +# context = AlgorithmContextCallback( +# dualGap = options.dualGap, +# κs = options.κs, +# shift = ψ.xk + ψ.sj, +# s_k_unshifted = s_k_unshifted, +# mk = ModelFunction(similar(x0), ψ), +# ) +# return iR2Solver( +# xk, +# ∇fk, +# mν∇fk, +# ψ, +# xkn, +# s, +# has_bnds, +# l_bound, +# u_bound, +# l_bound_m_x, +# u_bound_m_x, +# Fobj_hist, +# Hobj_hist, +# Complex_hist, +# callback_pointer, +# context, +# ) +# end + +function iR2Solver(reg_nlp::AbstractRegularizedNLPModel{T, V}; max_iter::Int = 10000) where {T, V} + x0 = reg_nlp.model.meta.x0 + l_bound = reg_nlp.model.meta.lvar + u_bound = reg_nlp.model.meta.uvar + + xk = similar(x0) + ∇fk = similar(x0) + mν∇fk = similar(x0) + xkn = similar(x0) + s = zero(x0) + has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf)) + if has_bnds + l_bound_m_x = similar(xk) + u_bound_m_x = similar(xk) + @. l_bound_m_x = l_bound - x0 + @. u_bound_m_x = u_bound - x0 + else + l_bound_m_x = similar(xk, 0) + u_bound_m_x = similar(xk, 0) + end + Fobj_hist = zeros(T, max_iter + 2) + Hobj_hist = zeros(T, max_iter + 2) + Complex_hist = zeros(Int, max_iter + 2) + + ψ = shifted(reg_nlp.h, xk) + # has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : + # shifted(reg_nlp.h, xk) + + return iR2Solver( + xk, + ∇fk, + mν∇fk, + ψ, + xkn, + s, + has_bnds, + l_bound, + u_bound, + l_bound_m_x, + u_bound_m_x, + Fobj_hist, + Hobj_hist, + Complex_hist, + ) +end + +""" + iR2(reg_nlp; kwargs…) + +A first-order quadratic regularization method for the problem + + min f(x) + h(x) + +where f: ℝⁿ → ℝ has a Lipschitz-continuous gradient, and h: ℝⁿ → ℝ is +lower semi-continuous, proper and prox-bounded. + +About each iterate xₖ, a step sₖ is computed as a solution of + + min φ(s; xₖ) + ½ σₖ ‖s‖² + ψ(s; xₖ) + +where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation of f about xₖ, +ψ(s; xₖ) is either h(xₖ + s) or an approximation of h(xₖ + s), ‖⋅‖ is a user-defined norm and σₖ > 0 is the regularization parameter. + +For advanced usage, first define a solver "R2Solver" to preallocate the memory used in the algorithm, and then call `solve!`: + + solver = R2Solver(reg_nlp) + solve!(solver, reg_nlp) + + stats = RegularizedExecutionStats(reg_nlp) + solver = R2Solver(reg_nlp) + solve!(solver, reg_nlp, stats) + +# Arguments +* `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the problem to solve, see `RegularizedProblems.jl`, `NLPModels.jl`. + +# Keyword arguments +- `x::V = nlp.meta.x0`: the initial guess; +- `atol::T = √eps(T)`: absolute tolerance; +- `rtol::T = √eps(T)`: relative tolerance; +- `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance +- `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); +- `max_time::Float64 = 30.0`: maximum time limit in seconds; +- `max_iter::Int = 10000`: maximum number of iterations; +- `verbose::Int = 0`: if > 0, display iteration details every `verbose` iteration; +- `σmin::T = eps(T)`: minimum value of the regularization parameter; +- `η1::T = √√eps(T)`: very successful iteration threshold; +- `η2::T = T(0.9)`: successful iteration threshold; +- `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ; +- `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful. + +The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure. + +# Output +The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. + +# Callback +The callback is called at each iteration. +The expected signature of the callback is `callback(nlp, solver, stats)`, and its output is ignored. +Changing any of the input arguments will affect the subsequent iterations. +In particular, setting `stats.status = :user` will stop the algorithm. +All relevant information should be available in `nlp` and `solver`. +Notably, you can access, and modify, the following: +- `solver.xk`: current iterate; +- `solver.∇fk`: current gradient; +- `stats`: structure holding the output of the algorithm (`GenericExecutionStats`), which contains, among other things: + - `stats.iter`: current iteration counter; + - `stats.objective`: current objective function value; + - `stats.solver_specific[:smooth_obj]`: current value of the smooth part of the objective function + - `stats.solver_specific[:nonsmooth_obj]`: current value of the nonsmooth part of the objective function + - `stats.status`: current status of the algorithm. Should be `:unknown` unless the algorithm has attained a stopping criterion. Changing this to anything will stop the algorithm, but you should use `:user` to properly indicate the intention. + - `stats.elapsed_time`: elapsed time in seconds. +""" +function iR2( + nlp::AbstractNLPModel{R, V}, + h, + options::ROSolverOptions{R}; + kwargs..., +) where {R <: Real, V} + kwargs_dict = Dict(kwargs...) + selected = pop!(kwargs_dict, :selected, 1:(nlp.meta.nvar)) + x0 = pop!(kwargs_dict, :x0, nlp.meta.x0) + reg_nlp = RegularizedNLPModel(nlp, h, selected) + return iR2( + reg_nlp, + x = x0, + atol = options.ϵa, + rtol = options.ϵr, + neg_tol = options.neg_tol, + verbose = options.verbose, + max_iter = options.maxIter, + max_time = options.maxTime, + σmin = options.σmin, + η1 = options.η1, + η2 = options.η2, + ν = options.ν, + γ = options.γ, + ) +end + +function iR2( + f::F, + ∇f!::G, + h::H, + options::ROSolverOptions{R}, + x0::AbstractVector{R}; + selected::AbstractVector{<:Integer} = 1:length(x0), + kwargs..., +) where {F <: Function, G <: Function, H, R <: Real} + nlp = FirstOrderModel(f, ∇f!, x0) + reg_nlp = RegularizedNLPModel(nlp, h, selected) + stats = iR2( + reg_nlp, + x = x0, + atol = options.ϵa, + rtol = options.ϵr, + neg_tol = options.neg_tol, + verbose = options.verbose, + max_iter = options.maxIter, + max_time = options.maxTime, + σmin = options.σmin, + η1 = options.η1, + η2 = options.η2, + ν = options.ν, + γ = options.γ, + ) + outdict = Dict( + :Fhist => stats.solver_specific[:Fhist], + :Hhist => stats.solver_specific[:Hhist], + :Chist => stats.solver_specific[:SubsolverCounter], + :ItersProx => stats.solver_specific[:ItersProx], + :NonSmooth => h, + :status => stats.status, + :fk => stats.solver_specific[:smooth_obj], + :hk => stats.solver_specific[:nonsmooth_obj], + :ξ => stats.solver_specific[:xi], + :elapsed_time => stats.elapsed_time, + ) + return stats.solution, stats.iter, outdict +end + +function iR2( + f::F, + ∇f!::G, + h::H, + options::ROSolverOptions{R}, + x0::AbstractVector{R}, + l_bound::AbstractVector{R}, + u_bound::AbstractVector{R}; + selected::AbstractVector{<:Integer} = 1:length(x0), + kwargs..., +) where {F <: Function, G <: Function, H, R <: Real} + nlp = FirstOrderModel(f, ∇f!, x0, lcon = l_bound, ucon = u_bound) + reg_nlp = RegularizedNLPModel(nlp, h, selected) + stats = iR2( + reg_nlp, + x = x0, + atol = options.ϵa, + rtol = options.ϵr, + neg_tol = options.neg_tol, + verbose = options.verbose, + max_iter = options.maxIter, + max_time = options.maxTime, + σmin = options.σmin, + η1 = options.η1, + η2 = options.η2, + ν = options.ν, + γ = options.γ, + ) + outdict = Dict( + :Fhist => stats.solver_specific[:Fhist], + :Hhist => stats.solver_specific[:Hhist], + :Chist => stats.solver_specific[:SubsolverCounter], + :ItersProx => stats.solver_specific[:ItersProx], + :NonSmooth => h, + :status => stats.status, + :fk => stats.solver_specific[:smooth_obj], + :hk => stats.solver_specific[:nonsmooth_obj], + :ξ => stats.solver_specific[:xi], + :elapsed_time => stats.elapsed_time, + ) + return stats.solution, stats.iter, outdict +end + +function iR2(reg_nlp::AbstractRegularizedNLPModel; kwargs...) + + # if h has exact prox, switch to R2 + # if !(shifted(reg_nlp.h, reg_nlp.model.meta.x0) isa InexactShiftedProximableFunction) + # @warn "h has exact prox, switching to R2" + # return R2(reg_nlp; kwargs...) + # end + + kwargs_dict = Dict(kwargs...) + max_iter = pop!(kwargs_dict, :max_iter, 10000) + + solver = iR2Solver(reg_nlp, max_iter = max_iter) + stats = GenericExecutionStats(reg_nlp.model) + cb = + (nlp, solver, stats) -> begin + solver.Fobj_hist[stats.iter + 1] = stats.solver_specific[:smooth_obj] + solver.Hobj_hist[stats.iter + 1] = stats.solver_specific[:nonsmooth_obj] + solver.Complex_hist[stats.iter + 1] += 1 + end + + solve!(solver, reg_nlp, stats; callback = cb, max_iter = max_iter, kwargs_dict...) + set_solver_specific!(stats, :Fhist, solver.Fobj_hist[1:(stats.iter + 1)]) + set_solver_specific!(stats, :Hhist, solver.Hobj_hist[1:(stats.iter + 1)]) + set_solver_specific!(stats, :SubsolverCounter, solver.Complex_hist[1:(stats.iter + 1)]) + return stats +end + +function SolverCore.solve!( + solver::iR2Solver{T}, + reg_nlp::AbstractRegularizedNLPModel{T, V}, + stats::GenericExecutionStats{T, V}; + callback = (args...) -> nothing, + x::V = reg_nlp.model.meta.x0, + atol::T = √eps(T), + rtol::T = √eps(T), + neg_tol::T = eps(T)^(1 / 4), + verbose::Int = 0, + max_iter::Int = 10000, + max_time::Float64 = 30.0, + max_eval::Int = -1, + σmin::T = eps(T), + η1::T = √√eps(T), + η2::T = T(0.9), + ν::T = eps(T)^(1 / 5), + γ::T = T(3), +) where {T, V} + reset!(stats) + + # Retrieve workspace + selected = reg_nlp.selected + h = reg_nlp.h + nlp = reg_nlp.model + + xk = solver.xk .= x + + # Make sure ψ has the correct shift + shift!(solver.ψ, xk) + + # ∇fk = solver.∇fk + mν∇fk = solver.mν∇fk + ψ = solver.ψ + κs = one(T) + dualGap = zero(T) + if ψ isa InexactShiftedProximableFunction + let ctx = ψ.h.context + κs = ctx.κs + dualGap = ctx.dualGap + end + end + + xkn = solver.xkn + s = solver.s + has_bnds = solver.has_bnds + if has_bnds + l_bound = solver.l_bound + u_bound = solver.u_bound + l_bound_m_x = solver.l_bound_m_x + u_bound_m_x = solver.u_bound_m_x + end + + # initialize parameters + improper = false + hk = @views h(xk[selected]) + if hk == Inf + verbose > 0 && @info "iR2: finding initial guess where nonsmooth term is finite" + prox!(xk, h, xk, one(eltype(xk))) + hk = @views h(xk[selected]) + hk < Inf || error("prox computation must be erroneous") + verbose > 0 && @debug "R2: found point where h has value" hk + end + improper = (hk == -Inf) + + if verbose > 0 + @info log_header( + [:iter, :fx, :hx, :norm_s_νInv, :ρ, :σ, :normx, :norms, :dualgap, :arrow], + [Int, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Char], + hdr_override = Dict{Symbol, String}( # TODO: Add this as constant dict elsewhere + :iter => "iter", + :fx => "f(x)", + :hx => "h(x)", + :norm_s_νInv => "‖s‖/ν", + :ρ => "ρ", + :σ => "σ", + :normx => "‖x‖", + :norms => "‖s‖", + :objgap => "dualGap_iR2", + :arrow => " ", + ), + colsep = 1, + ) + end + + local ξ::T + local ρk::T + σk = max(1 / ν, σmin) + ν = 1 / σk + sqrt_ξ_νInv = one(T) + norm_s_νInv = one(T) + + fk = obj(nlp, xk) + grad!(nlp, xk, solver.∇fk) + @. mν∇fk = -ν * solver.∇fk + + set_iter!(stats, 0) + start_time = time() + set_time!(stats, 0.0) + set_objective!(stats, fk + hk) + set_solver_specific!(stats, :smooth_obj, fk) + set_solver_specific!(stats, :nonsmooth_obj, hk) + + φk(d) = dot(solver.∇fk, d) + mk(d)::T = φk(d) + ψ(d)::T + + update_prox_context!(solver, stats, ψ) + prox!(s, ψ, mν∇fk, ν) + + mks = mk(s) + + # Record inner-level prox call counters (exclude feasibility prox above) + if ψ isa InexactShiftedProximableFunction + ctx = ψ.h.context + ctx.prox_calls_inner += 1 + ctx.prox_iters_inner += ctx.last_prox_iters + end + + ξ = hk - mks + max(1, abs(hk)) * 10 * eps() + + sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) + norm_s_νInv = norm(s) / ν + atol += rtol * norm_s_νInv # make stopping test absolute and relative + + # solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol * (√κs)) + solved = norm_s_νInv ≤ atol * κs # new stopping test according to theory + (ξ < 0 && sqrt_ξ_νInv > neg_tol) && + @warn "iR2: first prox-gradient step should produce a decrease but ξ = $(ξ) and sqrt_ξ_νInv = $(sqrt_ξ_νInv) > $(neg_tol)" + + set_solver_specific!(stats, :norm_s_νInv, norm_s_νInv) + set_status!( + stats, + get_status( + reg_nlp, + elapsed_time = stats.elapsed_time, + iter = stats.iter, + optimal = solved, + improper = improper, + max_eval = max_eval, + max_time = max_time, + max_iter = max_iter, + ), + ) + + callback(nlp, solver, stats) + + done = stats.status != :unknown + if done + println("iR2: terminating with ‖s‖/ν = $(norm_s_νInv) at iteration $(stats.iter)") + end + + while !done + + # Update xk, sigma_k + xkn .= xk .+ s + fkn = obj(nlp, xkn) + hkn = @views h(xkn[selected]) + improper = (hkn == -Inf) + + Δobj = (fk + hk) - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps() + ρk = Δobj / ξ + + verbose > 0 && + stats.iter % verbose == 0 && + @info log_row( + Any[ + stats.iter, + fk, + hk, + norm_s_νInv, + ρk, + σk, + norm(xk), + norm(s), + dualGap, + (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "="), + ], + colsep = 1, + ) + + if η1 ≤ ρk < Inf + xk .= xkn + if has_bnds + @. l_bound_m_x = l_bound - xk + @. u_bound_m_x = u_bound - xk + set_bounds!(ψ, l_bound_m_x, u_bound_m_x) + end + fk = fkn + hk = hkn + grad!(nlp, xk, solver.∇fk) + shift!(ψ, xk) + end + + if η2 ≤ ρk < Inf + σk = max(σk / γ, σmin) + end + if ρk < η1 || ρk == Inf + σk = σk * γ + end + + ν = 1 / σk + @. mν∇fk = -ν * solver.∇fk + + set_objective!(stats, fk + hk) + set_solver_specific!(stats, :smooth_obj, fk) + set_solver_specific!(stats, :nonsmooth_obj, hk) + set_iter!(stats, stats.iter + 1) + set_time!(stats, time() - start_time) + + # prepare callback context and pointer to callback function + update_prox_context!(solver, stats, ψ) + prox!(s, ψ, mν∇fk, ν) + mks = mk(s) + + # Record inner-level prox call counters per iteration + if ψ isa InexactShiftedProximableFunction + ctx = ψ.h.context + ctx.prox_calls_inner += 1 + ctx.prox_iters_inner += ctx.last_prox_iters + end + + ξ = hk - mks + max(1, abs(hk)) * 10 * eps() + + sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) + norm_s_νInv = norm(s) / ν + # solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol * (√κs)) + solved = norm_s_νInv ≤ atol * κs # new stopping test according to theory + (ξ < 0 && sqrt_ξ_νInv > neg_tol) && + @warn "iR2: prox-gradient step should produce a decrease but ξ = $(ξ) and sqrt_ξ_νInv = $(sqrt_ξ_νInv) > $(neg_tol)" + + set_solver_specific!(stats, :norm_s_νInv, norm_s_νInv) + set_status!( + stats, + get_status( + reg_nlp, + elapsed_time = stats.elapsed_time, + iter = stats.iter, + optimal = solved, + improper = improper, + max_eval = max_eval, + max_time = max_time, + max_iter = max_iter, + ), + ) + + callback(nlp, solver, stats) + + done = stats.status != :unknown + end + if verbose > 0 && stats.status == :first_order + @info log_row( + Any[ + stats.iter, + fk, + hk, + norm_s_νInv, + ρk, + σk, + norm(xk), + norm(s), + dualGap, + (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "="), + ], + colsep = 1, + ) + @info "iR2: terminating with ‖s‖/ν = $(norm_s_νInv)" + end + + if ψ isa InexactShiftedProximableFunction + set_solver_specific!(stats, :ItersProx, ψ.h.context.prox_stats[3]) + # Detailed inner prox stats (per-call means) + let ctx = ψ.h.context + total_calls_inner = ctx.prox_calls_inner + total_iters_inner = ctx.prox_iters_inner + mean_inner = total_calls_inner > 0 ? total_iters_inner / total_calls_inner : zero(eltype(xk)) + set_solver_specific!(stats, :total_iters_prox_inner, total_iters_inner) + set_solver_specific!(stats, :prox_calls_inner, total_calls_inner) + set_solver_specific!(stats, :mean_iters_prox_inner, mean_inner) + end + end + set_solution!(stats, xk) + set_residuals!(stats, zero(eltype(xk)), norm_s_νInv) + return stats +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 071616c0..306d0d68 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,7 +8,8 @@ using ADNLPModels, NLPModelsModifiers, RegularizedProblems, RegularizedOptimization, - SolverCore + SolverCore, + ProxTV Random.seed!(0) const global compound = 1 @@ -17,6 +18,7 @@ const global options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e const global bpdn, bpdn_nls, sol = bpdn_model(compound) const global bpdn2, bpdn_nls2, sol2 = bpdn_model(compound, bounds = true) const global λ = norm(grad(bpdn, zeros(bpdn.meta.nvar)), Inf) / 10 +const global p_val = 1.2 # for TVp regularization in iR2/iR2N include("test_AL.jl") @@ -134,5 +136,25 @@ for (mod, mod_name) ∈ ( end end +# iR2/iR2N +context = ProxTVContext(bpdn.meta.nvar, :lp, p_val, κs = 0.9, λ = 0.1) +hp = NormLp(0.1, p_val, context) +@testset "bpdn-iR2N-Lp" begin + x0 = zeros(bpdn.meta.nvar) + out = iR2N(LBFGSModel(bpdn), hp, options, x0 = x0) + @test typeof(out.solution) == typeof(bpdn.meta.x0) + @test length(out.solution) == bpdn.meta.nvar + @test typeof(out.dual_feas) == eltype(out.solution) + @test out.status == :first_order +end +@testset "bpdn-iR2-Lp" begin + x0 = zeros(bpdn.meta.nvar) + out = iR2(bpdn, hp, options, x0 = x0) + @test typeof(out.solution) == typeof(bpdn.meta.x0) + @test length(out.solution) == bpdn.meta.nvar + @test typeof(out.dual_feas) == eltype(out.solution) + @test out.status == :first_order +end + include("test_bounds.jl") include("test_allocs.jl") diff --git a/test/test_allocs.jl b/test/test_allocs.jl index 591fb58f..9ee92be9 100644 --- a/test/test_allocs.jl +++ b/test/test_allocs.jl @@ -60,6 +60,15 @@ end @test stats.status == :first_order end end + @testset "iR2" begin + reg_nlp = RegularizedNLPModel(bpdn, h) + solver = iR2Solver(reg_nlp) + stats = RegularizedExecutionStats(reg_nlp) + @test @wrappedallocs( + solve!(solver, reg_nlp, stats, ν = 1.0, atol = 1e-6, rtol = 1e-6, verbose = 0) + ) == 0 + @test stats.status == :first_order + end @testset "Augmented Lagrangian" begin continue # FIXME : fails due to type instabilities in ADNLPModels... reg_nlp = RegularizedNLPModel(hs8(backend = :generic), h)