simplify TR

Author: MaxenceGollier

src/TR_alg.jl

@@ -18,11 +18,6 @@ mutable struct TRSolver{
   s::V
   v0::V
   v1::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
@@ -36,39 +31,20 @@ function TRSolver(
   m_monotone::Int = 1,
 ) where {T, V, X}
   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)
-  has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf))
-  if has_bnds || subsolver == TRDHSolver
-    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)
 
   v0 = [(-1.0)^i for i = 0:(reg_nlp.model.meta.nvar - 1)]
   v0 ./= sqrt(reg_nlp.model.meta.nvar)
   v1 = similar(v0)
 
-  ψ =
-    has_bnds || subsolver == TRDHSolver ?
-    shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) :
-    shifted(reg_nlp.h, xk, T(1), χ)
-
-  Bk = hess_op(reg_nlp, xk)
-  sub_nlp = QuadraticModel(∇fk, Bk, x0 = x0)
   subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, χ = χ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
@@ -83,11 +59,6 @@ function TRSolver(
     s,
     v0,
     v1,
-    has_bnds,
-    l_bound,
-    u_bound,
-    l_bound_m_x,
-    u_bound_m_x,
     m_fh_hist,
     subsolver,
     subpb,
@@ -247,16 +218,12 @@ function SolverCore.solve!(
   s = solver.s
   χ = solver.χ
   m_fh_hist = solver.m_fh_hist .= T(-Inf)
-  has_bnds = solver.has_bnds
+  has_bnds = has_bounds(nlp)
 
   m_monotone = length(m_fh_hist) + 1
 
   set_radius!(solver.subpb, Δk)
   if has_bnds || isa(solver.subsolver, TRDHSolver) #TODO elsewhere ?
-    #l_bound_m_x, u_bound_m_x = solver.l_bound_m_x, solver.u_bound_m_x
-    #l_bound, u_bound = solver.l_bound, solver.u_bound
-    #update_bounds!(l_bound_m_x, u_bound_m_x, false, l_bound, u_bound, xk, Δk)
-    #set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
     set_bounds!(solver.subsolver.ψ, ψ.l, ψ.u)
   else
     set_radius!(solver.subsolver.ψ, Δk)
@@ -328,20 +295,8 @@ function SolverCore.solve!(
   set_solver_specific!(stats, :prox_evals, prox_evals + 1)
   m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)
 
-  # models
-  φ1 = let ∇fk = ∇fk
-    d -> dot(∇fk, d)
-  end
-  mk1 = let ψ = ψ, φ1 = φ1
-    d -> φ1(d) + ψ(d)
-  end
-
-  mk = let ψ = ψ, solver = solver
-    d -> obj(solver.subpb.model, d) + ψ(d)::T
-  end
-
   prox!(s, ψ, mν∇fk, ν₁)
-  ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
+  ξ1 = hk - obj(mk, s, cauchy = true) + max(1, abs(hk)) * 10 * eps()
   ξ1 > 0 || error("TR: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)")
   sqrt_ξ1_νInv = sqrt(ξ1 / ν₁)
 
@@ -372,14 +327,13 @@ function SolverCore.solve!(
     sub_atol = stats.iter == 0 ? 1e-5 : max(sub_atol, min(1e-2, sqrt_ξ1_νInv))
     ∆_effective = min(β * χ(s), Δk)
 
-    if has_bnds || isa(solver.subsolver, TRDHSolver) #TODO elsewhere ?
-      update_bounds!(l_bound_m_x, u_bound_m_x, false, l_bound, u_bound, xk, Δk)
-      set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
-      set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
+    set_radius!(solver.subpb, ∆_effective)
+    if has_bnds || isa(solver.subsolver, TRDHSolver)
+      set_bounds!(solver.subsolver.ψ, ψ.l, ψ.u)
     else
       set_radius!(solver.subsolver.ψ, ∆_effective)
-      set_radius!(ψ, ∆_effective)
     end
+
     with_logger(subsolver_logger) do
       if isa(solver.subsolver, TRDHSolver) #FIXME
         solver.subsolver.D.d[1] = 1/ν₁
@@ -415,7 +369,7 @@ function SolverCore.solve!(
 
     fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk
     Δobj = fhmax - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
-    ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()
+    ξ = hk - obj(mk, s) + max(1, abs(hk)) * 10 * eps()
 
     if (ξ ≤ 0 || isnan(ξ))
       error("TR: failed to compute a step: ξ = $ξ")
@@ -444,25 +398,20 @@ function SolverCore.solve!(
 
     if η2 ≤ ρk < Inf
       Δk = max(Δk, γ * sNorm)
+      set_radius!(solver.subpb, Δk)
       if !(has_bnds || isa(solver.subsolver, TRDHSolver))
-        set_radius!(ψ, Δk)
         set_radius!(solver.subsolver.ψ, Δk)
       end
     end
 
     if η1 ≤ ρk < Inf
       xk .= xkn
-      if has_bnds || isa(solver.subsolver, TRDHSolver)
-        update_bounds!(l_bound_m_x, u_bound_m_x, false, l_bound, u_bound, xk, Δk)
-        set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
-        set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
-      end
+      
+      shift!(solver.subpb, xk, compute_grad = compute_grad)
+
       fk = fkn
       hk = hkn
 
-      shift!(ψ, xk)
-      grad!(nlp, xk, ∇fk)
-
       if quasiNewtTest
         @. ∇fk⁻ = ∇fk - ∇fk⁻
         push!(nlp, s, ∇fk⁻) # update QN operator
@@ -481,12 +430,10 @@ function SolverCore.solve!(
 
     if ρk < η1 || ρk == Inf
       Δk = Δk / 2
+      set_radius!(solver.subpb, Δk)
       if has_bnds || isa(solver.subsolver, TRDHSolver)
-        update_bounds!(l_bound_m_x, u_bound_m_x, false, l_bound, u_bound, xk, Δk)
-        set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
-        set_bounds!(solver.subsolver.ψ, l_bound_m_x, u_bound_m_x)
+        set_bounds!(solver.subsolver.ψ, ψ.l, ψ.u)
       else
-        set_radius!(ψ, Δk)
         set_radius!(solver.subsolver.ψ, Δk)
       end
       set_step_status!(stats, :rejected)
@@ -505,7 +452,7 @@ function SolverCore.solve!(
     @. mν∇fk = -ν₁ * ∇fk
 
     prox!(s, ψ, mν∇fk, ν₁)
-    ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps()
+    ξ1 = hk - obj(mk, s, cauchy = true) + max(1, abs(hk)) * 10 * eps()
     sqrt_ξ1_νInv = sqrt(ξ1 / ν₁)
 
     solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol)