Introduction of R2N Solver in JSOSolvers

src/R2N.jl

@@ -7,6 +7,7 @@ struct ShiftedLBFGSSolver <: AbstractShiftedLBFGSSolver
   # Shifted LBFGS-specific fields
 end
 
+
 """
     R2N(nlp; kwargs...)
 
@@ -147,8 +148,9 @@ function SolverCore.reset!(solver::R2NSolver{T}) where {T}
 end
 function SolverCore.reset!(solver::R2NSolver{T}, nlp::AbstractNLPModel) where {T}
   fill!(solver.obj_vec, typemin(T))
-  @assert (length(solver.gn) == 0) || isa(nlp, QuasiNewtonModel)
+  # @assert (length(solver.gn) == 0) || isa(nlp, QuasiNewtonModel)
   solver.H = isa(nlp, QuasiNewtonModel) ? nlp.op : hess_op!(nlp, solver.x, solver.Hs)
+
   solver
 end
 
@@ -282,14 +284,23 @@ function SolverCore.solve!(
     ∇fk .*= -1
     subsolve!(solver, s, zero(T), n, subsolver_verbose)
 
-    slope = dot(s, ∇fk) # = -dot(s, ∇fk) but ∇fk is negative
+    slope = dot(s, ∇fk)
     mul!(Hs, H, s)
     curv = dot(s, Hs)
 
-    ΔTk = slope + curv / 2
+    ΔTk = slope - curv / 2
     ck .= x .+ s
     fck = obj(nlp, ck)
-
+    ϵ = eps(T)
+
+    if ΔTk <= 0
+      ΔTk += max(one(T), fck) * 10 * ϵ
+      if ΔTk <= 0
+        stats.status = :neg_pred
+        done = true
+        continue
+      end
+    end
     if non_mono_size > 1  #non-monotone behaviour
       k = mod(stats.iter, non_mono_size) + 1
       solver.obj_vec[k] = stats.objective
@@ -316,6 +327,7 @@ function SolverCore.solve!(
       norm_∇fk = norm(∇fk)
     else
       μk = μk * λ
+      ∇fk .*= -1
     end
 
     set_iter!(stats, stats.iter + 1)
@@ -386,12 +398,15 @@ function subsolve!(R2N::R2NSolver, s, atol, n, subsolver_verbose)
       H,
       ∇f, #b 
       λ = σ,
-      itmax = 2 * n,
+      itmax = max(2 * n, 50),
       atol = atol,
       rtol = cgtol,
       verbose = subsolver_verbose,
     )
     s .= subsolver_type.x
+    if norm(subsolver_type.x) < atol
+      println("X_norm is:" ,norm(subsolver_type.x)," the status is:"     ,subsolver_type.stats.status)
+    end
   elseif subsolver_type isa KrylovSolver
     Krylov.solve!(
       subsolver_type,
@@ -403,7 +418,6 @@ function subsolve!(R2N::R2NSolver, s, atol, n, subsolver_verbose)
       verbose = subsolver_verbose,
     )
     s .= subsolver_type.x
-
   elseif subsolver_type isa ShiftedLBFGSSolver
     solve_shifted_system!(s, H, ∇f, σ)
   else

test/allocs.jl

@@ -31,12 +31,14 @@ end
 if Sys.isunix()
   @testset "Allocation tests" begin
     @testset "$symsolver" for symsolver in
-                              (:LBFGSSolver, :FoSolver, :FomoSolver, :TrunkSolver, :TronSolver, :R2NSolver)
+                              (:LBFGSSolver, :FoSolver, :FomoSolver, :TrunkSolver, :TronSolver, :R2NSolver, :R2NExactSolver)
       for model in NLPModelsTest.nlp_problems
         nlp = eval(Meta.parse(model))()
         if unconstrained(nlp) || (bound_constrained(nlp) && (symsolver == :TronSolver))
           if (symsolver == :FoSolver || symsolver == :FomoSolver)
             solver = eval(symsolver)(nlp; M = 2) # nonmonotone configuration allocates extra memory
+          elseif symsolver == :R2NExactSolver
+            solver = eval(:R2NExactSolver)(LBFGSModel(nlp), subsolver_type = JSOSolvers.ShiftedLBFGSSolver)
           else
             solver = eval(symsolver)(nlp)
           end

test/consistency.jl

@@ -10,28 +10,44 @@ function consistency()
   @testset "Consistency" begin
     args = Pair{Symbol, Number}[:atol => 1e-6, :rtol => 1e-6, :max_eval => 20000, :max_time => 60.0]
 
-    @testset "NLP with $mtd" for mtd in [trunk, lbfgs, tron, R2, fomo, R2N]
+    @testset "NLP with $mtd" for (mtd, solver) in [
+            ("trunk", trunk),
+            ("lbfgs", lbfgs),
+            ("tron", tron),
+            ("R2", R2),
+            # ("R2N", R2N),
+            ("R2N_exact", (nlp; kwargs...) -> R2N(LBFGSModel(nlp), subsolver_type = JSOSolvers.ShiftedLBFGSSolver; kwargs...)),
+            ("fomo", fomo),
+          ]
       with_logger(NullLogger()) do
         reset!(unlp)
-        stats = mtd(unlp; args...)
+        stats = solver(unlp; args...)
         @test stats isa GenericExecutionStats
         @test stats.status == :first_order
         reset!(unlp)
-        stats = mtd(unlp; max_eval = 1)
+        stats = solver(unlp; max_eval = 1)
         @test stats.status == :max_eval
         slow_nlp = ADNLPModel(x -> begin
           sleep(0.1)
           f(x)
         end, unlp.meta.x0)
-        stats = mtd(slow_nlp; max_time = 0.0)
+        stats = solver(slow_nlp; max_time = 0.0)
         @test stats.status == :max_time
       end
     end
 
-    @testset "Quasi-Newton NLP with $mtd" for mtd in [trunk, lbfgs, tron, R2, fomo, R2N]
+    @testset "Quasi-Newton NLP with $mtd" for (mtd, solver) in [
+      ("trunk", trunk),
+      ("lbfgs", lbfgs),
+      ("tron", tron),
+      ("R2", R2),
+      # ("R2N", R2N),
+      ("R2N_exact", (nlp; kwargs...) -> R2N(LBFGSModel(nlp), subsolver_type = JSOSolvers.ShiftedLBFGSSolver; kwargs...)),
+      ("fomo", fomo),
+    ]
       with_logger(NullLogger()) do
         reset!(qnlp)
-        stats = mtd(qnlp; args...)
+        stats = solver(qnlp; args...)
         @test stats isa GenericExecutionStats
         @test stats.status == :first_order
       end

test/restart.jl

@@ -1,5 +1,6 @@
 @testset "Test restart with a different initial guess: $fun" for (fun, s) in (
-  (:R2N, :R2NSolver),
+  # (:R2N, :R2NSolver),
+  (:R2N_exact, :R2NSolver),
   (:R2, :FoSolver),
   (:fomo, :FomoSolver),
   (:lbfgs, :LBFGSSolver),
@@ -8,9 +9,15 @@
 )
   f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2
   nlp = ADNLPModel(f, [-1.2; 1.0])
+  if fun == :R2N_exact
+    nlp = LBFGSModel(nlp)
+    solver = eval(s)(nlp,subsolver_type = JSOSolvers.ShiftedLBFGSSolver)
+  else 
+    solver = eval(s)(nlp)
+  end
 
   stats = GenericExecutionStats(nlp)
-  solver = eval(s)(nlp)
+
   stats = SolverCore.solve!(solver, nlp, stats)
   @test stats.status == :first_order
   @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)
@@ -45,7 +52,8 @@ end
 end
 
 @testset "Test restart with a different problem: $fun" for (fun, s) in (
-  (:R2N, :R2NSolver),
+  # (:R2N, :R2NSolver),
+  (:R2N_exact, :R2NSolver),
   (:R2, :FoSolver),
   (:fomo, :FomoSolver),
   (:lbfgs, :LBFGSSolver),
@@ -54,15 +62,25 @@ end
 )
   f(x) = (x[1] - 1)^2 + 4 * (x[2] - x[1]^2)^2
   nlp = ADNLPModel(f, [-1.2; 1.0])
+  if fun == :R2N_exact
+    nlp = LBFGSModel(nlp)
+    solver = eval(s)(nlp,subsolver_type = JSOSolvers.ShiftedLBFGSSolver)
+  else 
+    solver = eval(s)(nlp) 
+  end
 
   stats = GenericExecutionStats(nlp)
-  solver = eval(s)(nlp)
   stats = SolverCore.solve!(solver, nlp, stats)
   @test stats.status == :first_order
   @test isapprox(stats.solution, [1.0; 1.0], atol = 1e-6)
 
   f2(x) = (x[1])^2 + 4 * (x[2] - x[1]^2)^2
   nlp = ADNLPModel(f2, [-1.2; 1.0])
+  if fun == :R2N_exact
+    nlp = LBFGSModel(nlp)
+  else 
+    solver = eval(s)(nlp) 
+  end
   SolverCore.reset!(solver, nlp)
 
   stats = SolverCore.solve!(solver, nlp, stats, atol = 1e-10, rtol = 1e-10)

test/test_solvers.jl

@@ -8,7 +8,8 @@ function tests()
         ("lbfgs", lbfgs),
         ("tron", tron),
         ("R2", R2),
-        ("R2N", R2N),
+        # ("R2N", R2N),
+        ("R2N_exact", (nlp; kwargs...) -> R2N(LBFGSModel(nlp), subsolver_type = JSOSolvers.ShiftedLBFGSSolver; kwargs...)),
         ("fomo_r2", fomo),
         ("fomo_tr", (nlp; kwargs...) -> fomo(nlp, step_backend = JSOSolvers.tr_step(); kwargs...)),
       ]