From f3fd5ef98e00271013f81be411a72e4290126780 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Tue, 24 Mar 2026 18:12:38 -0400
Subject: [PATCH 01/12] start subproblem construction update

---
 src/R2N.jl                     | 17 +++++++----------
 src/RegularizedOptimization.jl |  2 ++
 2 files changed, 9 insertions(+), 10 deletions(-)

diff --git a/src/R2N.jl b/src/R2N.jl
index acb05d94..f662bd24 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -33,6 +33,7 @@ end
 
 function R2NSolver(
   reg_nlp::AbstractRegularizedNLPModel{T, V};
+  subproblem = quadratic_subproblem,
   subsolver = R2Solver,
   m_monotone::Int = 1,
 ) where {T, V}
@@ -65,17 +66,13 @@ function R2NSolver(
   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)
+  subproblem = subproblem(reg_nlp, xk)
+  subsolver = subsolver(subproblem)
+  substats = RegularizedExecutionStats(subproblem)
 
-  Bk = hess_op(reg_nlp, xk)
-  sub_nlp = QuadraticModel(∇fk, Bk, σ = T(1), x0 = x0)
-  subpb = RegularizedNLPModel(sub_nlp, ψ)
-  substats = RegularizedExecutionStats(subpb)
-  subsolver = subsolver(subpb)
+  ψ = subproblem.h
 
-  return R2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subpb)}(
+  return R2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subproblem)}(
     xk,
     ∇fk,
     ∇fk⁻,
@@ -94,7 +91,7 @@ function R2NSolver(
     u_bound_m_x,
     m_fh_hist,
     subsolver,
-    subpb,
+    subproblem,
     substats,
   )
 end
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index de226752..3f645bcb 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -19,6 +19,8 @@ using Percival: AugLagModel, update_y!, update_μ!
 
 import SolverCore.reset!
 
+abstract type AbstractRegularizedOptimizationSolver <: AbstractOptimizationSolver end
+
 const callback_docstring = "
 The callback is called at each iteration.
 The expected signature of the callback is `callback(reg_nlp, solver, stats)`, and its output is ignored.

From f3a77b7422ddd009291cb1af08dbb7efa634bcba Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Tue, 24 Mar 2026 18:22:38 -0400
Subject: [PATCH 02/12] add subproblem constructors

---
 src/utils.jl | 16 ++++++++++++++++
 1 file changed, 16 insertions(+)

diff --git a/src/utils.jl b/src/utils.jl
index 15c53d9e..880809df 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -189,3 +189,19 @@ function update_bounds!(
     @. u_bound_m_x = min(Δ, u_bound - xk)
   end
 end
+
+# Given min_x f(x) + h(x), return min_s φ(s; x) + ψ(s; x), where φ(s; x) = f(x) + ∇f(x)ᵀs + ½ sᵀBs
+function quadratic_subproblem(
+  reg_nlp::AbstractRegularizedNLPModel{T, V},
+  x::V,
+) where{T, V}
+  #TODO
+end
+
+# Given min_x f(x) + h(x), return min_s φ(s; x) + ψ(s; x), where φ(s; x) = f(x) + ∇f(x)ᵀs
+function linear_subproblem(
+  reg_nlp::AbstractRegularizedNLPModel{T, V},
+  x::V,
+) where{T, V}
+  #TODO
+end
\ No newline at end of file

From a2929b9399903e8db19cdc3abf943570e33da4da Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 25 Mar 2026 10:27:04 -0400
Subject: [PATCH 03/12] add a subproblem type

---
 src/R2N.jl                     | 47 +++++---------------------------
 src/RegularizedOptimization.jl |  2 ++
 src/subproblem.jl              | 50 ++++++++++++++++++++++++++++++++++
 src/utils.jl                   | 16 -----------
 4 files changed, 59 insertions(+), 56 deletions(-)
 create mode 100644 src/subproblem.jl

diff --git a/src/R2N.jl b/src/R2N.jl
index f662bd24..b0219282 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -20,11 +20,6 @@ mutable struct R2NSolver{
   s1::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
@@ -33,13 +28,10 @@ end
 
 function R2NSolver(
   reg_nlp::AbstractRegularizedNLPModel{T, V};
-  subproblem = quadratic_subproblem,
   subsolver = R2Solver,
   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)
@@ -54,19 +46,12 @@ function R2NSolver(
   v0 ./= sqrt(reg_nlp.model.meta.nvar)
   v1 = similar(v0)
 
-  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
+  l_bound_m_x = has_bounds(reg_nlp.model) ? similar(x0) : nothing
+  u_bound_m_x = has_bounds(reg_nlp.model) ? similar(x0) : nothing
+
   m_fh_hist = fill(T(-Inf), m_monotone - 1)
 
-  subproblem = subproblem(reg_nlp, xk)
+  subproblem = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, l_bound_m_x = l_bound_m_x, u_bound_m_x = u_bound_m_x, ∇f = ∇fk)
   subsolver = subsolver(subproblem)
   substats = RegularizedExecutionStats(subproblem)
 
@@ -84,11 +69,6 @@ function R2NSolver(
     s1,
     v0,
     v1,
-    has_bnds,
-    l_bound,
-    u_bound,
-    l_bound_m_x,
-    u_bound_m_x,
     m_fh_hist,
     subsolver,
     subproblem,
@@ -242,7 +222,7 @@ function SolverCore.solve!(
   xk = solver.xk .= x
 
   # Make sure ψ has the correct shift 
-  shift!(solver.ψ, xk)
+  shift!(solver.subpb, xk, compute_grad = compute_grad)
 
   ∇fk = solver.∇fk
   ∇fk⁻ = solver.∇fk⁻
@@ -252,14 +232,7 @@ function SolverCore.solve!(
   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, u_bound = solver.l_bound, solver.u_bound
-    l_bound_m_x, u_bound_m_x = solver.l_bound_m_x, solver.u_bound_m_x
-    update_bounds!(l_bound_m_x, u_bound_m_x, l_bound, u_bound, xk)
-    set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
-  end
   m_monotone = length(m_fh_hist) + 1
 
   # initialize parameters
@@ -298,7 +271,6 @@ function SolverCore.solve!(
   local prox_evals::Int = 0
 
   fk = compute_obj ? obj(nlp, xk) : stats.solver_specific[:smooth_obj]
-  compute_grad && grad!(nlp, xk, ∇fk)
   qn_copy!(nlp, solver, stats)
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
@@ -446,17 +418,12 @@ function SolverCore.solve!(
 
     if η1 ≤ ρk < Inf
       xk .= xkn
-      if has_bnds
-        update_bounds!(l_bound_m_x, u_bound_m_x, l_bound, u_bound, xk)
-        set_bounds!(ψ, l_bound_m_x, u_bound_m_x)
-      end
+      shift!(solver.subpb, xk)
+
       #update functions
       fk = fkn
       hk = hkn
 
-      shift!(ψ, xk)
-      grad!(nlp, xk, ∇fk)
-
       if quasiNewtTest
         qn_update_y!(nlp, solver, stats)
         push!(nlp, s, solver.y)
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index 3f645bcb..ab88da40 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -17,6 +17,7 @@ using LinearOperators,
   SolverCore
 using Percival: AugLagModel, update_y!, update_μ!
 
+import ShiftedProximalOperators.shift!
 import SolverCore.reset!
 
 abstract type AbstractRegularizedOptimizationSolver <: AbstractOptimizationSolver end
@@ -41,6 +42,7 @@ Notably, you can access, and modify, the following:
 
 include("utils.jl")
 include("input_struct.jl")
+include("subproblem.jl")
 include("TR_alg.jl")
 include("TRDH_alg.jl")
 include("R2_alg.jl")
diff --git a/src/subproblem.jl b/src/subproblem.jl
new file mode 100644
index 00000000..115bc120
--- /dev/null
+++ b/src/subproblem.jl
@@ -0,0 +1,50 @@
+abstract type AbstractShiftedProximableNLPModel{T, V} <: AbstractRegularizedNLPModel{T, V} end
+
+struct ShiftedProximableQuadraticNLPModel{T, V, M <: AbstractNLPModel{T, V}, H <: ShiftedProximableFunction, I, P <: AbstractRegularizedNLPModel{T, V}} <:
+       AbstractShiftedProximableNLPModel{T, V}
+  model::M
+  h::H
+  selected::I
+  parent::P
+end
+
+function ShiftedProximableQuadraticNLPModel(
+  reg_nlp::AbstractRegularizedNLPModel{T, V}, 
+  x::V;
+  l_bound_m_x::VN = nothing,
+  u_bound_m_x::VN = nothing,
+  ∇f::VNG = nothing,
+) where {T, V, VN <: Union{V, Nothing}, VNG <: Union{V, Nothing}}
+  nlp, h, selected = reg_nlp.model, reg_nlp.h, reg_nlp.selected
+
+  @assert !(has_bounds(nlp) && isnothing(l_bound_m_x) && isnothing(u_bound_m_x)) 
+    "RegularizedOptimization: bounds are required for the quadratic subproblem when the NLP has bounds."
+
+  # FIXME: `shifted` call ignores the `selected` argument when there are no bounds!
+  ψ = has_bounds(nlp) ? shifted(h, x, l_bound_m_x, u_bound_m_x, selected) : shifted(h, x)
+
+  B = hess_op(reg_nlp, x)
+  isnothing(∇f) && (∇f = grad(nlp, x))
+  φ = QuadraticModel(∇f, B, x0 = x, regularize = true)
+
+  ShiftedProximableQuadraticNLPModel(φ, ψ, selected, reg_nlp)
+end
+
+function shift!(
+  reg_nlp::ShiftedProximableQuadraticNLPModel{T, V},
+  x::V;
+  compute_grad::Bool = true
+) where{T, V}
+  nlp, h = reg_nlp.parent.model, reg_nlp.parent.h
+  φ, ψ = reg_nlp.model, reg_nlp.h
+
+  if has_bounds(nlp)
+    update_bounds!(ψ.l, ψ.u, nlp.meta.lvar, nlp.meta.uvar, x)
+  end
+  shift!(ψ, x)
+
+  g = φ.data.c
+  compute_grad && grad!(nlp, x, g)
+
+  # The hessian is implicitely updated since it was defined as hess_op(nlp, x)
+end
\ No newline at end of file
diff --git a/src/utils.jl b/src/utils.jl
index 880809df..b10422e8 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -188,20 +188,4 @@ function update_bounds!(
     @. l_bound_m_x = max(-Δ, l_bound - xk)
     @. u_bound_m_x = min(Δ, u_bound - xk)
   end
-end
-
-# Given min_x f(x) + h(x), return min_s φ(s; x) + ψ(s; x), where φ(s; x) = f(x) + ∇f(x)ᵀs + ½ sᵀBs
-function quadratic_subproblem(
-  reg_nlp::AbstractRegularizedNLPModel{T, V},
-  x::V,
-) where{T, V}
-  #TODO
-end
-
-# Given min_x f(x) + h(x), return min_s φ(s; x) + ψ(s; x), where φ(s; x) = f(x) + ∇f(x)ᵀs
-function linear_subproblem(
-  reg_nlp::AbstractRegularizedNLPModel{T, V},
-  x::V,
-) where{T, V}
-  #TODO
 end
\ No newline at end of file

From 02b684ae876cb06e8c4ea73de9fe828d404ea624 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 25 Mar 2026 10:35:15 -0400
Subject: [PATCH 04/12] fix diff

---
 src/R2N.jl                     | 12 ++++++------
 src/RegularizedOptimization.jl |  2 --
 2 files changed, 6 insertions(+), 8 deletions(-)

diff --git a/src/R2N.jl b/src/R2N.jl
index b0219282..ccf4833c 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -51,13 +51,13 @@ function R2NSolver(
 
   m_fh_hist = fill(T(-Inf), m_monotone - 1)
 
-  subproblem = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, l_bound_m_x = l_bound_m_x, u_bound_m_x = u_bound_m_x, ∇f = ∇fk)
-  subsolver = subsolver(subproblem)
-  substats = RegularizedExecutionStats(subproblem)
+  subpb= ShiftedProximableQuadraticNLPModel(reg_nlp, xk, l_bound_m_x = l_bound_m_x, u_bound_m_x = u_bound_m_x, ∇f = ∇fk)
+  subsolver = subsolver(subpb)
+  substats = RegularizedExecutionStats(subpb)
 
-  ψ = subproblem.h
+  ψ = subpb.h
 
-  return R2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subproblem)}(
+  return R2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subpb)}(
     xk,
     ∇fk,
     ∇fk⁻,
@@ -71,7 +71,7 @@ function R2NSolver(
     v1,
     m_fh_hist,
     subsolver,
-    subproblem,
+    subpb,
     substats,
   )
 end
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index ab88da40..895f782e 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -20,8 +20,6 @@ using Percival: AugLagModel, update_y!, update_μ!
 import ShiftedProximalOperators.shift!
 import SolverCore.reset!
 
-abstract type AbstractRegularizedOptimizationSolver <: AbstractOptimizationSolver end
-
 const callback_docstring = "
 The callback is called at each iteration.
 The expected signature of the callback is `callback(reg_nlp, solver, stats)`, and its output is ignored.

From 56163b2f9409f315049a41e1673fbcf233fd2613 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 25 Mar 2026 10:56:21 -0400
Subject: [PATCH 05/12] add a doc

---
 src/subproblem.jl | 54 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 53 insertions(+), 1 deletion(-)

diff --git a/src/subproblem.jl b/src/subproblem.jl
index 115bc120..27219c17 100644
--- a/src/subproblem.jl
+++ b/src/subproblem.jl
@@ -1,5 +1,38 @@
 abstract type AbstractShiftedProximableNLPModel{T, V} <: AbstractRegularizedNLPModel{T, V} end
 
+"""
+    subproblem = ShiftedProximableQuadraticNLPModel(reg_nlp, x; kwargs...)
+
+Given a regularized NLP model `reg_nlp` representing the problem
+
+    minimize f(x) + h(x),
+  
+construct a shifted quadratic model around `x`:
+
+    minimize  φ(s; x) + ½ σ ‖s‖² + ψ(s; x),
+
+where φ(s ; x) = f(x) + ∇f(x)ᵀs + ½ sᵀBs 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.
+
+The ShiftedProximableQuadraticNLPModel is made of the following components:
+
+- `model <: AbstractNLPModel`: represents φ, the quadratic approximation of the smooth part of the objective function;
+- `h <: ShiftedProximableFunction`: represents ψ, the shifted version of the nonsmooth part of the model;
+- `selected`: the subset of variables to which the regularizer h should be applied (default: all).
+- `parent`: the original regularized NLP model from which the subproblem was derived.
+
+# Arguments
+- `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the regularized NLP model for which the subproblem is being constructed.
+- `x::V`: the point around which the quadratic model is constructed.
+
+# Keyword Arguments
+- `l_bound_m_x::VN = nothing`: the vector of lower bounds minus `x` (i.e., l - x), required if the original NLP model has bounds.
+- `u_bound_m_x::VN = nothing`: the vector of upper bounds minus `x` (i.e., u - x), required if the original NLP model has bounds.
+- `∇f::VNG = nothing`: the gradient of the smooth part of the objective function at `x`. If not provided, it will be computed.
+
+The matrix B is constructed as a `LinearOperator` and is the returned value of `hess_op(reg_nlp, x)` (see https://jso.dev/NLPModels.jl/stable/reference/#NLPModels.hess_op`).
+φ is constructed as a `QuadraticModel`, (see https://github.com/JuliaSmoothOptimizers/QuadraticModels.jl).
+"""
 struct ShiftedProximableQuadraticNLPModel{T, V, M <: AbstractNLPModel{T, V}, H <: ShiftedProximableFunction, I, P <: AbstractRegularizedNLPModel{T, V}} <:
        AbstractShiftedProximableNLPModel{T, V}
   model::M
@@ -30,6 +63,25 @@ function ShiftedProximableQuadraticNLPModel(
   ShiftedProximableQuadraticNLPModel(φ, ψ, selected, reg_nlp)
 end
 
+"""
+    shift!(reg_nlp::ShiftedProximableQuadraticNLPModel, x; compute_grad = true)
+
+Update the shifted quadratic model `reg_nlp` at the point `x`. 
+i.e. given the shifted quadratic model around `y`:
+
+    minimize  φ(s; y) + ½ σ ‖s‖² + ψ(s; y),
+
+update it to be around `x`:
+
+    minimize  φ(s; x) + ½ σ ‖s‖² + ψ(s; x).
+
+# Arguments
+- `reg_nlp::ShiftedProximableQuadraticNLPModel`: the shifted quadratic model to be updated.
+- `x::V`: the point around which the shifted quadratic model should be updated.
+
+# Keyword Arguments
+- `compute_grad::Bool = true`: whether the gradient of the smooth part of the model should be updated.
+"""
 function shift!(
   reg_nlp::ShiftedProximableQuadraticNLPModel{T, V},
   x::V;
@@ -46,5 +98,5 @@ function shift!(
   g = φ.data.c
   compute_grad && grad!(nlp, x, g)
 
-  # The hessian is implicitely updated since it was defined as hess_op(nlp, x)
+  # The hessian is implicitly updated since it was defined as hess_op(nlp, x)
 end
\ No newline at end of file

From 6fb91b2399bc1e287e47b1f26c547fb7b1f99970 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 25 Mar 2026 16:02:14 -0400
Subject: [PATCH 06/12] simplify R2N

---
 src/R2N.jl                     |  32 +++--------
 src/RegularizedOptimization.jl |   2 -
 src/subproblem.jl              | 102 ---------------------------------
 3 files changed, 7 insertions(+), 129 deletions(-)
 delete mode 100644 src/subproblem.jl

diff --git a/src/R2N.jl b/src/R2N.jl
index ccf4833c..c16a967c 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -46,12 +46,9 @@ function R2NSolver(
   v0 ./= sqrt(reg_nlp.model.meta.nvar)
   v1 = similar(v0)
 
-  l_bound_m_x = has_bounds(reg_nlp.model) ? similar(x0) : nothing
-  u_bound_m_x = has_bounds(reg_nlp.model) ? similar(x0) : nothing
-
   m_fh_hist = fill(T(-Inf), m_monotone - 1)
 
-  subpb= ShiftedProximableQuadraticNLPModel(reg_nlp, xk, l_bound_m_x = l_bound_m_x, u_bound_m_x = u_bound_m_x, ∇f = ∇fk)
+  subpb= ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk)
   subsolver = subsolver(subpb)
   substats = RegularizedExecutionStats(subpb)
 
@@ -301,26 +298,10 @@ 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)
 
-  φ1 = let ∇fk = ∇fk
-    d -> dot(∇fk, d)
-  end
-
-  mk1 = let ψ = ψ
-    d -> φ1(d) + ψ(d)::T
-  end
-
-  mk = let ψ = ψ, model = solver.subpb.model
-    d -> begin
-      temp_σ = model.data.σ
-      model.data.σ = zero(T)
-      _obj = obj(model, d) + ψ(d)
-      model.data.σ = temp_σ
-      return _obj
-    end
-  end
+  mk = solver.subpb
 
   prox!(s1, ψ, mν∇fk, ν₁)
-  mks = mk1(s1)
+  mks = obj(mk, s1, cauchy = true, skip_sigma = true)
 
   ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps()
   sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁)
@@ -350,7 +331,7 @@ function SolverCore.solve!(
   while !done
     sub_atol = stats.iter == 0 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5), sqrt_ξ1_νInv * 1e-3)
 
-    solver.subpb.model.data.σ = σk
+    update_sigma!(solver.subpb, σk)
     isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1/ν₁)
     if isa(solver.subsolver, R2Solver) #FIXME
       solve!(
@@ -384,7 +365,7 @@ function SolverCore.solve!(
     xkn .= xk .+ s
     fkn = obj(nlp, xkn)
     hkn = @views h(xkn[selected])
-    mks = mk(s)
+    mks = obj(mk, s, skip_sigma = true) 
 
     fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk
     Δobj = fhmax - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
@@ -418,6 +399,7 @@ function SolverCore.solve!(
 
     if η1 ≤ ρk < Inf
       xk .= xkn
+
       shift!(solver.subpb, xk)
 
       #update functions
@@ -462,7 +444,7 @@ function SolverCore.solve!(
 
     @. mν∇fk = - ν₁ * ∇fk
     prox!(s1, ψ, mν∇fk, ν₁)
-    mks = mk1(s1)
+    mks = obj(mk, s1, cauchy = true,skip_sigma = true)
 
     ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps()
 
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index 895f782e..de226752 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -17,7 +17,6 @@ using LinearOperators,
   SolverCore
 using Percival: AugLagModel, update_y!, update_μ!
 
-import ShiftedProximalOperators.shift!
 import SolverCore.reset!
 
 const callback_docstring = "
@@ -40,7 +39,6 @@ Notably, you can access, and modify, the following:
 
 include("utils.jl")
 include("input_struct.jl")
-include("subproblem.jl")
 include("TR_alg.jl")
 include("TRDH_alg.jl")
 include("R2_alg.jl")
diff --git a/src/subproblem.jl b/src/subproblem.jl
deleted file mode 100644
index 27219c17..00000000
--- a/src/subproblem.jl
+++ /dev/null
@@ -1,102 +0,0 @@
-abstract type AbstractShiftedProximableNLPModel{T, V} <: AbstractRegularizedNLPModel{T, V} end
-
-"""
-    subproblem = ShiftedProximableQuadraticNLPModel(reg_nlp, x; kwargs...)
-
-Given a regularized NLP model `reg_nlp` representing the problem
-
-    minimize f(x) + h(x),
-  
-construct a shifted quadratic model around `x`:
-
-    minimize  φ(s; x) + ½ σ ‖s‖² + ψ(s; x),
-
-where φ(s ; x) = f(x) + ∇f(x)ᵀs + ½ sᵀBs 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.
-
-The ShiftedProximableQuadraticNLPModel is made of the following components:
-
-- `model <: AbstractNLPModel`: represents φ, the quadratic approximation of the smooth part of the objective function;
-- `h <: ShiftedProximableFunction`: represents ψ, the shifted version of the nonsmooth part of the model;
-- `selected`: the subset of variables to which the regularizer h should be applied (default: all).
-- `parent`: the original regularized NLP model from which the subproblem was derived.
-
-# Arguments
-- `reg_nlp::AbstractRegularizedNLPModel{T, V}`: the regularized NLP model for which the subproblem is being constructed.
-- `x::V`: the point around which the quadratic model is constructed.
-
-# Keyword Arguments
-- `l_bound_m_x::VN = nothing`: the vector of lower bounds minus `x` (i.e., l - x), required if the original NLP model has bounds.
-- `u_bound_m_x::VN = nothing`: the vector of upper bounds minus `x` (i.e., u - x), required if the original NLP model has bounds.
-- `∇f::VNG = nothing`: the gradient of the smooth part of the objective function at `x`. If not provided, it will be computed.
-
-The matrix B is constructed as a `LinearOperator` and is the returned value of `hess_op(reg_nlp, x)` (see https://jso.dev/NLPModels.jl/stable/reference/#NLPModels.hess_op`).
-φ is constructed as a `QuadraticModel`, (see https://github.com/JuliaSmoothOptimizers/QuadraticModels.jl).
-"""
-struct ShiftedProximableQuadraticNLPModel{T, V, M <: AbstractNLPModel{T, V}, H <: ShiftedProximableFunction, I, P <: AbstractRegularizedNLPModel{T, V}} <:
-       AbstractShiftedProximableNLPModel{T, V}
-  model::M
-  h::H
-  selected::I
-  parent::P
-end
-
-function ShiftedProximableQuadraticNLPModel(
-  reg_nlp::AbstractRegularizedNLPModel{T, V}, 
-  x::V;
-  l_bound_m_x::VN = nothing,
-  u_bound_m_x::VN = nothing,
-  ∇f::VNG = nothing,
-) where {T, V, VN <: Union{V, Nothing}, VNG <: Union{V, Nothing}}
-  nlp, h, selected = reg_nlp.model, reg_nlp.h, reg_nlp.selected
-
-  @assert !(has_bounds(nlp) && isnothing(l_bound_m_x) && isnothing(u_bound_m_x)) 
-    "RegularizedOptimization: bounds are required for the quadratic subproblem when the NLP has bounds."
-
-  # FIXME: `shifted` call ignores the `selected` argument when there are no bounds!
-  ψ = has_bounds(nlp) ? shifted(h, x, l_bound_m_x, u_bound_m_x, selected) : shifted(h, x)
-
-  B = hess_op(reg_nlp, x)
-  isnothing(∇f) && (∇f = grad(nlp, x))
-  φ = QuadraticModel(∇f, B, x0 = x, regularize = true)
-
-  ShiftedProximableQuadraticNLPModel(φ, ψ, selected, reg_nlp)
-end
-
-"""
-    shift!(reg_nlp::ShiftedProximableQuadraticNLPModel, x; compute_grad = true)
-
-Update the shifted quadratic model `reg_nlp` at the point `x`. 
-i.e. given the shifted quadratic model around `y`:
-
-    minimize  φ(s; y) + ½ σ ‖s‖² + ψ(s; y),
-
-update it to be around `x`:
-
-    minimize  φ(s; x) + ½ σ ‖s‖² + ψ(s; x).
-
-# Arguments
-- `reg_nlp::ShiftedProximableQuadraticNLPModel`: the shifted quadratic model to be updated.
-- `x::V`: the point around which the shifted quadratic model should be updated.
-
-# Keyword Arguments
-- `compute_grad::Bool = true`: whether the gradient of the smooth part of the model should be updated.
-"""
-function shift!(
-  reg_nlp::ShiftedProximableQuadraticNLPModel{T, V},
-  x::V;
-  compute_grad::Bool = true
-) where{T, V}
-  nlp, h = reg_nlp.parent.model, reg_nlp.parent.h
-  φ, ψ = reg_nlp.model, reg_nlp.h
-
-  if has_bounds(nlp)
-    update_bounds!(ψ.l, ψ.u, nlp.meta.lvar, nlp.meta.uvar, x)
-  end
-  shift!(ψ, x)
-
-  g = φ.data.c
-  compute_grad && grad!(nlp, x, g)
-
-  # The hessian is implicitly updated since it was defined as hess_op(nlp, x)
-end
\ No newline at end of file

From 529d9c4d053705e9119a963c8a6e7e6845cac179 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 25 Mar 2026 16:28:28 -0400
Subject: [PATCH 07/12] further simplify R2N

---
 src/R2N.jl | 36 +++++++++++++++---------------------
 1 file changed, 15 insertions(+), 21 deletions(-)

diff --git a/src/R2N.jl b/src/R2N.jl
index c16a967c..125b02e3 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -4,7 +4,6 @@ import SolverCore.solve!
 
 mutable struct R2NSolver{
   T <: Real,
-  G <: ShiftedProximableFunction,
   V <: AbstractVector{T},
   ST <: AbstractOptimizationSolver,
   PB <: AbstractRegularizedNLPModel,
@@ -14,7 +13,6 @@ mutable struct R2NSolver{
   ∇fk⁻::V
   y::V
   mν∇fk::V
-  ψ::G
   xkn::V
   s::V
   s1::V
@@ -52,15 +50,12 @@ function R2NSolver(
   subsolver = subsolver(subpb)
   substats = RegularizedExecutionStats(subpb)
 
-  ψ = subpb.h
-
-  return R2NSolver{T, typeof(ψ), V, typeof(subsolver), typeof(subpb)}(
+  return R2NSolver{T, V, typeof(subsolver), typeof(subpb)}(
     xk,
     ∇fk,
     ∇fk⁻,
     y,
     mν∇fk,
-    ψ,
     xkn,
     s,
     s1,
@@ -183,7 +178,7 @@ function R2N(reg_nlp::AbstractRegularizedNLPModel; kwargs...)
 end
 
 function SolverCore.solve!(
-  solver::R2NSolver{T, G, V},
+  solver::R2NSolver{T, V},
   reg_nlp::AbstractRegularizedNLPModel{T, V},
   stats::GenericExecutionStats{T, V};
   callback = (args...) -> nothing,
@@ -208,7 +203,7 @@ function SolverCore.solve!(
   compute_obj::Bool = true,
   compute_grad::Bool = true,
   sub_kwargs::NamedTuple = NamedTuple(),
-) where {T, V, G}
+) where {T, V}
   reset!(stats)
 
   # Retrieve workspace
@@ -218,13 +213,14 @@ function SolverCore.solve!(
 
   xk = solver.xk .= x
 
-  # Make sure ψ has the correct shift 
-  shift!(solver.subpb, xk, compute_grad = compute_grad)
+  mk = solver.subpb
+  φ, ψ = mk.model, mk.h
+  
+  shift!(mk, xk, compute_grad = compute_grad)
 
   ∇fk = solver.∇fk
   ∇fk⁻ = solver.∇fk⁻
   mν∇fk = solver.mν∇fk
-  ψ = solver.ψ
   xkn = solver.xkn
   s = solver.s
   s1 = solver.s1
@@ -275,9 +271,9 @@ function SolverCore.solve!(
   found_λ = true
 
   if opnorm_maxiter ≤ 0
-    λmax, found_λ = opnorm(solver.subpb.model.data.H)
+    λmax, found_λ = opnorm(φ.data.H)
   else
-    λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
+    λmax = power_method!(φ.data.H, solver.v0, solver.v1, opnorm_maxiter)
   end
   found_λ || error("operator norm computation failed")
 
@@ -298,8 +294,6 @@ 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)
 
-  mk = solver.subpb
-
   prox!(s1, ψ, mν∇fk, ν₁)
   mks = obj(mk, s1, cauchy = true, skip_sigma = true)
 
@@ -413,9 +407,9 @@ function SolverCore.solve!(
       end
 
       if opnorm_maxiter ≤ 0
-        λmax, found_λ = opnorm(solver.subpb.model.data.H)
+        λmax, found_λ = opnorm(φ.data.H)
       else
-        λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
+        λmax = power_method!(φ.data.H, solver.v0, solver.v1, opnorm_maxiter)
       end
       found_λ || error("operator norm computation failed")
       set_step_status!(stats, :accepted)
@@ -487,16 +481,16 @@ end
 
 function _qn_grad_update_y!(
   nlp::AbstractNLPModel{T, V},
-  solver::R2NSolver{T, G, V},
+  solver::R2NSolver{T, V},
   stats::GenericExecutionStats,
-) where {T, V, G}
+) where {T, V}
   @. solver.y = solver.∇fk - solver.∇fk⁻
 end
 
 function _qn_grad_copy!(
   nlp::AbstractNLPModel{T, V},
-  solver::R2NSolver{T, G, V},
+  solver::R2NSolver{T, V},
   stats::GenericExecutionStats,
-) where {T, V, G}
+) where {T, V}
   solver.∇fk⁻ .= solver.∇fk
 end

From 10b14eead4248cf5186231ee80f783e575b919fa Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Mon, 30 Mar 2026 13:47:59 -0400
Subject: [PATCH 08/12] start TR simplification

---
 src/R2N.jl    |  2 +-
 src/TR_alg.jl | 28 +++++++++++++---------------
 2 files changed, 14 insertions(+), 16 deletions(-)

diff --git a/src/R2N.jl b/src/R2N.jl
index 125b02e3..3633e89b 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -46,7 +46,7 @@ function R2NSolver(
 
   m_fh_hist = fill(T(-Inf), m_monotone - 1)
 
-  subpb= ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk)
+  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk)
   subsolver = subsolver(subpb)
   substats = RegularizedExecutionStats(subpb)
 
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 7ebec003..8916d7a0 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -4,7 +4,6 @@ import SolverCore.solve!
 
 mutable struct TRSolver{
   T <: Real,
-  G <: ShiftedProximableFunction,
   V <: AbstractVector{T},
   N,
   ST <: AbstractOptimizationSolver,
@@ -14,7 +13,6 @@ mutable struct TRSolver{
   ∇fk::V
   ∇fk⁻::V
   mν∇fk::V
-  ψ::G
   χ::N
   xkn::V
   s::V
@@ -71,7 +69,7 @@ function TRSolver(
 
   Bk = hess_op(reg_nlp, xk)
   sub_nlp = QuadraticModel(∇fk, Bk, x0 = x0)
-  subpb = RegularizedNLPModel(sub_nlp, ψ)
+  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, χ = χ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
 
@@ -80,7 +78,6 @@ function TRSolver(
     ∇fk,
     ∇fk⁻,
     mν∇fk,
-    ψ,
     χ,
     xkn,
     s,
@@ -206,7 +203,7 @@ function TR(reg_nlp::AbstractRegularizedNLPModel{T, V}; kwargs...) where {T, V}
 end
 
 function SolverCore.solve!(
-  solver::TRSolver{T, G, V},
+  solver::TRSolver{T, V},
   reg_nlp::AbstractRegularizedNLPModel{T, V},
   stats::GenericExecutionStats{T, V};
   callback = (args...) -> nothing,
@@ -228,7 +225,7 @@ function SolverCore.solve!(
   opnorm_maxiter::Int = 5,
   compute_obj::Bool = true,
   compute_grad::Bool = true,
-) where {T, G, V}
+) where {T, V}
   reset!(stats)
 
   # Retrieve workspace
@@ -238,13 +235,14 @@ function SolverCore.solve!(
 
   xk = solver.xk .= x
 
-  # Make sure ψ has the correct shift 
-  shift!(solver.ψ, xk)
+  mk = solver.subpb
+  φ, ψ = mk.model, mk.h
+  
+  shift!(mk, xk, compute_grad = compute_grad)
 
   ∇fk = solver.∇fk
   ∇fk⁻ = solver.∇fk⁻
   mν∇fk = solver.mν∇fk
-  ψ = solver.ψ
   xkn = solver.xkn
   s = solver.s
   χ = solver.χ
@@ -253,14 +251,14 @@ function SolverCore.solve!(
 
   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_bound_m_x, u_bound_m_x)
+    #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!(ψ, Δk)
     set_radius!(solver.subsolver.ψ, Δk)
   end
 

From 982af7e4c18d49eba687dccec5982b7639d62cdf Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Mon, 30 Mar 2026 17:42:57 -0400
Subject: [PATCH 09/12] simplify TR

---
 src/TR_alg.jl | 81 +++++++++------------------------------------------
 1 file changed, 14 insertions(+), 67 deletions(-)

diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 8916d7a0..2de0d778 100644
--- a/src/TR_alg.jl
+++ b/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,8 +31,6 @@ 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)
@@ -45,16 +38,6 @@ function TRSolver(
   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)
 
@@ -62,13 +45,6 @@ function TRSolver(
   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)

From 77d24f9d29869a095fdbdbdac3ca2701629968b5 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Tue, 31 Mar 2026 09:11:44 -0400
Subject: [PATCH 10/12] fix TRDH subsolver

---
 src/TR_alg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 2de0d778..557690b3 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -45,7 +45,7 @@ function TRSolver(
   v0 ./= sqrt(reg_nlp.model.meta.nvar)
   v1 = similar(v0)
 
-  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, χ = χ)
+  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, χ = χ, has_bounds = (has_bounds(reg_nlp.model) || subsolver == TRDHSolver))
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
 

From a8c211207fbf31941945a88480b602e5411a0020 Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Tue, 31 Mar 2026 10:19:36 -0400
Subject: [PATCH 11/12] update subsolver constructor logic

---
 src/TR_alg.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 557690b3..82d9b064 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -45,7 +45,8 @@ function TRSolver(
   v0 ./= sqrt(reg_nlp.model.meta.nvar)
   v1 = similar(v0)
 
-  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, χ = χ, has_bounds = (has_bounds(reg_nlp.model) || subsolver == TRDHSolver))
+  indicator_type = has_bounds(reg_nlp.model) || (subsolver == TRDHSolver) ? :box : :ball
+  subpb = ShiftedProximableQuadraticNLPModel(reg_nlp, xk, ∇f = ∇fk, indicator_type = indicator_type, tr_norm = χ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
 
@@ -407,7 +408,7 @@ function SolverCore.solve!(
     if η1 ≤ ρk < Inf
       xk .= xkn
       
-      shift!(solver.subpb, xk, compute_grad = compute_grad)
+      shift!(solver.subpb, xk)
 
       fk = fkn
       hk = hkn

From 7a9bd99770fdd15e296801baf60ab90a91df1c7b Mon Sep 17 00:00:00 2001
From: MaxenceGollier <mgollier45@gmail.com>
Date: Wed, 1 Apr 2026 15:02:09 -0400
Subject: [PATCH 12/12] update name

---
 src/R2N.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/R2N.jl b/src/R2N.jl
index 3633e89b..785901f2 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -325,7 +325,7 @@ function SolverCore.solve!(
   while !done
     sub_atol = stats.iter == 0 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5), sqrt_ξ1_νInv * 1e-3)
 
-    update_sigma!(solver.subpb, σk)
+    set_sigma!(solver.subpb, σk)
     isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1/ν₁)
     if isa(solver.subsolver, R2Solver) #FIXME
       solve!(