add updated parameter of subsolvers and add sigam-k to input.jl structure

Author: MohamedLaghdafHABIBOULLAH

src/R2DH.jl

@@ -61,6 +61,8 @@ function R2DH(
       hess_op(nlp, x0),
       options,
       x0;
+      l_bound = nlp.meta.lvar,
+      u_bound = nlp.meta.uvar,
       kwargs...,
     )
     sqrt_ξ_νInv = outdict[:sqrt_ξ_νInv]
@@ -99,6 +101,7 @@ function R2DH(
   maxIter = options.maxIter
   maxTime = options.maxTime
   σmin = options.σmin
+  σk = options.σk
   η1 = options.η1
   η2 = options.η2
   ν = options.ν
@@ -148,6 +151,7 @@ function R2DH(
   time_hist = zeros(maxIter+1)
   FHobj_hist = fill!(Vector{R}(undef, Mmonotone), R(-Inf))
   Complex_hist = zeros(Int, maxIter+1)
+  k_prox = 0
   if verbose > 0
     #! format: off
     @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %1s" "iter" "f(x)" "h(x)" "√(ξ/ν)" "ρ" "σ" "‖x‖" "‖s‖" ""
@@ -156,13 +160,13 @@ function R2DH(
 
   local ξ
   k = 0
-  σk = σmin
 
   fk = f(xk)
   ∇fk = similar(xk)
   ∇f!(∇fk, xk)
   ∇fk⁻ = copy(∇fk) 
   spectral_test = isa(D, SpectralGradient)
+ # D.d .= D.d .+ σk
   Dkσk = D.d .+ σk
   DNorm = norm(D.d, Inf)
 
@@ -174,15 +178,7 @@ function R2DH(
   tired = maxIter > 0 && k ≥ maxIter || elapsed_time > maxTime
 
   while !(optimal || tired)
-    k = k + 1
-    elapsed_time = time() - start_time
-    Fobj_hist[k] = fk
-    Hobj_hist[k] = hk
-    time_hist[k] = elapsed_time
-    Mmonotone > 0 && (FHobj_hist[mod(k-1, Mmonotone) + 1] = fk + hk)
-
-    #Dkσk .= max.(Dkσk, eps(R)) 
-
+    k_prox += 1
     # model with diagonal hessian
     φ(d) = ∇fk' * d + (d' * (Dkσk .* d)) / 2
     mk(d) = φ(d) + ψ(d)
@@ -203,7 +199,15 @@ function R2DH(
       continue
     end
 
-    Complex_hist[k] += 1
+    k = k + 1
+    elapsed_time = time() - start_time
+    Fobj_hist[k] = fk
+    Hobj_hist[k] = hk
+    time_hist[k] = elapsed_time
+    Mmonotone > 0 && (FHobj_hist[mod(k-1, Mmonotone) + 1] = fk + hk)
+
+    Complex_hist[k] += k_prox
+    k_prox = 0
     xkn .= xk .+ s
     fkn = f(xkn)
     hkn = h(xkn[selected])
@@ -260,6 +264,7 @@ function R2DH(
     end
 
     Dkσk .= D.d .+ σk
+    DNorm = norm(D.d, Inf)
     ν = 1 / ((DNorm + σk) * (1 + θ))
 
     tired = maxIter > 0 && k ≥ maxIter

src/R2N.jl

@@ -51,6 +51,7 @@ function R2N(
   subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(),
   subsolver = R2,
   subsolver_options = ROSolverOptions(ϵa = options.ϵa),
+  Mmonotone::Int = 0, 
   selected::AbstractVector{<:Integer} = 1:(f.meta.nvar),
 ) where {H, R}
   start_time = time()
@@ -71,6 +72,7 @@ function R2N(
   σmin = options.σmin
   α = options.α
   β = options.β
+  σk = options.σk
 
   # store initial values of the subsolver_options fields that will be modified
   ν_subsolver = subsolver_options.ν
@@ -94,7 +96,6 @@ function R2N(
 
   # initialize parameters
   #σk = max(1 / options.ν, σmin) #SVM
-  σk = σmin
   xk = copy(x0)
   hk = h(xk[selected])
   if hk == Inf
@@ -112,6 +113,7 @@ function R2N(
 
   Fobj_hist = zeros(maxIter)
   Hobj_hist = zeros(maxIter)
+  FHobj_hist = fill!(Vector{R}(undef, Mmonotone), R(-Inf))
   Complex_hist = zeros(Int, maxIter)
   if verbose > 0
     #! format: off
@@ -143,6 +145,7 @@ function R2N(
     elapsed_time = time() - start_time
     Fobj_hist[k] = fk
     Hobj_hist[k] = hk
+    Mmonotone > 0 && (FHobj_hist[mod(k-1, Mmonotone) + 1] = fk + hk)
 
     # model for first prox-gradient step and ξ1
     φ1(d) = ∇fk' * d
@@ -168,9 +171,6 @@ function R2N(
     ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
     ξ1 > 0 || error("R2N: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
     sqrt_ξ1_νInv = sqrt(ξ1 * νInv)
-    # println("sqrt_ξ1_νInv: ", sqrt_ξ1_νInv)
-    # println("ξ1: ", ξ1)
-    # println("νInv: ", νInv)
 
     if ξ1 ≥ 0 && k == 1
       ϵ_increment = ϵr * sqrt_ξ1_νInv
@@ -186,9 +186,9 @@ function R2N(
     s1 = copy(s)
 
   #  subsolver_options.ϵa = k == 1 ? 1.0e-1 : max(ϵ_subsolver, min(1.0e-2, ξ1 / 10))
-    subsolver_options.ϵa = k == 1 ? 1.0e-3 : max(ϵ_subsolver, min(1e-3, sqrt_ξ1_νInv)) # 1.0e-5 default
+    subsolver_options.ϵa = k == 1 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5) , sqrt_ξ1_νInv * 1e-3) # 1.0e-5 default
     @debug "setting inner stopping tolerance to" subsolver_options.optTol
-    subsolver_args = subsolver == R2DH ? (SpectralGradient(1., f.meta.nvar),) : ()
+    subsolver_args = subsolver == R2DH ? (SpectralGradient(νInv, f.meta.nvar),) : ()
     s, iter, _ = with_logger(subsolver_logger) do
       subsolver(φ, ∇φ!, ψ, subsolver_args..., subsolver_options, s)
     end
@@ -208,18 +208,21 @@ function R2N(
     hkn = h(xkn[selected])
     hkn == -Inf && error("nonsmooth term is not proper")
     mks = mk(s) #- σk * dot(s, s) / 2
-    Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
+
+    fhmax = Mmonotone > 0 ? maximum(FHobj_hist) : fk + hk
+    Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps()
+    Δmod = fhmax - (fk + mks) + max(1, abs(hk)) * 10 * eps()
     ξ = hk - mks + max(1, abs(hk)) * 10 * eps()
 
     if (ξ ≤ 0 || isnan(ξ))
       error("R2N: failed to compute a step: ξ = $ξ")
     end
 
-    ρk = Δobj / ξ
+    ρk = Δobj / Δmod
 
     R2N_stat = (η2 ≤ ρk < Inf) ? "↗" : (ρk < η1 ? "↘" : "=")
 
-    if (verbose > 0) && (k % ptf == 0)
+    if (verbose > 0) && ((k % ptf == 0) || (k == 1))
       #! format: off
       @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt_ξ1_νInv sqrt(ξ1) ρk σk norm(xk) norm(s) λmax R2N_stat
       #! format: off
@@ -251,8 +254,7 @@ function R2N(
     if ρk < η1 || ρk == Inf
         σk = σk * γ
     end
-    νInv = (1 + θ) *( σk + λmax)
-
+    νInv = (1 + θ) * (σk + λmax)
     tired = k ≥ maxIter || elapsed_time > maxTime
   end
 

src/input_struct.jl

@@ -9,6 +9,7 @@ mutable struct ROSolverOptions{R}
   maxIter::Int  # maximum amount of inner iterations
   maxTime::Float64 #maximum time allotted to the algorithm in s
   σmin::R # minimum σk allowed for LM/R2 method
+  σk::R # initial σk
   η1::R  # step acceptance threshold
   η2::R  # trust-region increase threshold
   α::R  # νk Δ^{-1} parameter
@@ -27,6 +28,7 @@ mutable struct ROSolverOptions{R}
     maxIter::Int = 10000,
     maxTime::Float64 = 3600.0,
     σmin::R = eps(R),
+    σk::R = eps(R)^(1 / 5),
     η1::R = √√eps(R),
     η2::R = R(0.9),
     α::R = 1 / eps(R),
@@ -44,6 +46,7 @@ mutable struct ROSolverOptions{R}
     @assert maxIter ≥ 0
     @assert maxTime ≥ 0
     @assert σmin ≥ 0
+    @assert σk ≥ 0
     @assert 0 < η1 < η2 < 1
     @assert α > 0
     @assert ν > 0
@@ -59,6 +62,7 @@ mutable struct ROSolverOptions{R}
       maxIter,
       maxTime,
       σmin,
+      σk,
       η1,
       η2,
       α,