R2: remove nu keyword argument and replace with sigma

Author: MaxenceGollier

src/LMTR_alg.jl

@@ -197,7 +197,7 @@ function LMTR(
     set_bounds!(ψ, max.(-∆_effective, l_bound - xk), min.(∆_effective, u_bound - xk)) :
     set_radius!(ψ, ∆_effective)
     subsolver_options.Δk = ∆_effective / 10
-    subsolver_options.ν = ν
+    subsolver_options.σk = 1/ν
     subsolver_args = subsolver == TRDH ? (SpectralGradient(1 / ν, nls.meta.nvar),) : ()
     s, iter, _ = with_logger(subsolver_logger) do
       subsolver(φ, ∇φ!, ψ, subsolver_args..., subsolver_options, s)

src/LM_alg.jl

@@ -318,18 +318,15 @@ function SolverCore.solve!(
     sub_atol = stats.iter == 0 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5), sqrt_ξ1_νInv * 1e-3)
     solver.subpb.model.σ = σk
     isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1/ν)
-    if isa(solver.subsolver, R2Solver) #FIXME
-      solve!(solver.subsolver, solver.subpb, solver.substats, x = s, atol = sub_atol, ν = ν)
-    else
-      solve!(
-        solver.subsolver,
-        solver.subpb,
-        solver.substats,
-        x = s,
-        atol = sub_atol,
-        σk = σk, #FIXME
-      )
-    end
+
+    solve!(
+      solver.subsolver,
+      solver.subpb,
+      solver.substats,
+      x = s,
+      atol = sub_atol,
+      σk = 1/ν, #FIXME
+    )
 
     prox_evals += solver.substats.iter
     s .= solver.substats.solution

src/R2DH.jl

@@ -312,6 +312,7 @@ function SolverCore.solve!(
   ∇fk⁻ .= ∇fk
   spectral_test = isa(D, SpectralGradient)
 
+  σk = max(σk, σmin)
   @. dkσk = D.d .+ σk
   DNorm = norm(D.d, Inf)
 

src/R2N.jl

@@ -304,6 +304,7 @@ function SolverCore.solve!(
   λmax, found_λ = opnorm(solver.subpb.model.B)
   found_λ || error("operator norm computation failed")
 
+  σk = max(σk, σmin)
   ν₁ = θ / (λmax + σk)
 
   sqrt_ξ1_νInv = one(T)
@@ -366,27 +367,15 @@ function SolverCore.solve!(
 
     solver.subpb.model.σ = σk
     isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1/ν₁)
-    if isa(solver.subsolver, R2Solver) #FIXME
-      solve!(
-        solver.subsolver,
-        solver.subpb,
-        solver.substats;
-        x = s1,
-        ν = ν₁,
-        atol = sub_atol,
-        sub_kwargs...,
-      )
-    else
-      solve!(
-        solver.subsolver,
-        solver.subpb,
-        solver.substats;
-        x = s1,
-        σk = σk,
-        atol = sub_atol,
-        sub_kwargs...,
-      )
-    end
+    solve!(
+      solver.subsolver,
+      solver.subpb,
+      solver.substats;
+      x = s1,
+      σk = σk,
+      atol = sub_atol,
+      sub_kwargs...,
+    )
 
     prox_evals += solver.substats.iter
     s .= solver.substats.solution

src/R2_alg.jl

@@ -151,7 +151,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u
 - `σ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/σ;
+- `σk::T = eps(T)^(1 / 5)`: initial value of the regularization parameter;
 - `γ::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.
@@ -198,7 +198,7 @@ function R2(
     σmin = options.σmin,
     η1 = options.η1,
     η2 = options.η2,
-    ν = options.ν,
+    σk = options.σk,
     γ = options.γ;
     kwargs_dict...,
   )
@@ -237,7 +237,7 @@ function R2(
     σmin = options.σmin,
     η1 = options.η1,
     η2 = options.η2,
-    ν = options.ν,
+    σk = options.σk,
     γ = options.γ;
     kwargs...,
   )
@@ -280,7 +280,7 @@ function R2(
     σmin = options.σmin,
     η1 = options.η1,
     η2 = options.η2,
-    ν = options.ν,
+    σk = options.σk,
     γ = options.γ;
     kwargs...,
   )
@@ -335,7 +335,7 @@ function SolverCore.solve!(
   σmin::T = eps(T),
   η1::T = √√eps(T),
   η2::T = T(0.9),
-  ν::T = eps(T)^(1 / 5),
+  σk::T = eps(T)^(1 / 5),
   γ::T = T(3),
 ) where {T, V}
   reset!(stats)
@@ -396,7 +396,7 @@ function SolverCore.solve!(
 
   local ξ::T
   local ρk::T
-  σk = max(1 / ν, σmin)
+  σk = max(σk, σmin)
   ν = 1 / σk
   sqrt_ξ_νInv = one(T)
 

src/TR_alg.jl

@@ -366,7 +366,7 @@ function SolverCore.solve!(
           solver.substats;
           x = s,
           atol = stats.iter == 0 ? 1e-5 : max(sub_atol, min(1e-2, sqrt_ξ1_νInv)),
-          ν = ν₁,
+          σk = 1/ν₁,
           sub_kwargs...,
         )
       end

test/runtests.jl

@@ -12,7 +12,7 @@ using ADNLPModels,
 
 const global compound = 1
 const global nz = 10 * compound
-const global options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
+const global options = ROSolverOptions(σk = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
 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

test/test_allocs.jl

@@ -54,10 +54,7 @@ end
         reg_nlp = RegularizedNLPModel(LBFGSModel(bpdn), h)
         solver = eval(solver)(reg_nlp)
         stats = RegularizedExecutionStats(reg_nlp)
-        solver_name == "R2" &&
-          @test @wrappedallocs(solve!(solver, reg_nlp, stats, ν = 1.0, atol = 1e-6, rtol = 1e-6)) ==
-                0
-        solver_name == "R2DH" && @test @wrappedallocs(
+        (solver_name == "R2DH" || solver_name == "R2") && @test @wrappedallocs(
           solve!(solver, reg_nlp, stats, σk = 1.0, atol = 1e-6, rtol = 1e-6)
         ) == 0
         solver_name == "TRDH" &&