diff --git a/Project.toml b/Project.toml
index 82030991..00f4e5ed 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,6 +19,7 @@ NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
+QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
@@ -33,6 +34,7 @@ NLPModelsModifiers = "0.7"
 OptimizationProblems = "0.9.2"
 Percival = "0.7.2"
 ProximalOperators = "0.15"
+QuadraticModels = "0.9.15"
 RegularizedProblems = "0.1.3"
 ShiftedProximalOperators = "0.2"
 SolverCore = "0.3.10"
diff --git a/src/R2N.jl b/src/R2N.jl
index 51ad3c56..e9879109 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -19,6 +19,7 @@ mutable struct R2NSolver{
   s::V
   s1::V
   v0::V
+  v1::V
   has_bnds::Bool
   l_bound::V
   u_bound::V
@@ -50,6 +51,7 @@ function R2NSolver(
 
   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 = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf))
   if has_bnds
@@ -68,7 +70,7 @@ function R2NSolver(
     shifted(reg_nlp.h, xk)
 
   Bk = hess_op(reg_nlp, xk)
-  sub_nlp = R2NModel(Bk, ∇fk, T(1), x0)
+  sub_nlp = QuadraticModel(∇fk, Bk, σ = T(1), x0 = x0)
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
@@ -84,6 +86,7 @@ function R2NSolver(
     s,
     s1,
     v0,
+    v1,
     has_bnds,
     l_bound,
     u_bound,
@@ -98,7 +101,7 @@ end
 
 function SolverCore.reset!(solver::R2NSolver)
   _reset_power_method!(solver.v0)
-  B = solver.subpb.model.B
+  B = solver.subpb.model.data.H
   isa(B, AbstractLinearOperator) && LinearOperators.reset!(B)
 end
 
@@ -306,9 +309,9 @@ function SolverCore.solve!(
   found_λ = true
 
   if opnorm_maxiter ≤ 0
-    λmax, found_λ = opnorm(solver.subpb.model.B)
+    λmax, found_λ = opnorm(solver.subpb.model.data.H)
   else
-    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+    λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
   end
   found_λ || error("operator norm computation failed")
 
@@ -337,8 +340,14 @@ function SolverCore.solve!(
     d -> φ1(d) + ψ(d)::T
   end
 
-  mk = let ψ = ψ, solver = solver
-    d -> obj(solver.subpb.model, d, skip_sigma = true) + ψ(d)::T
+  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
 
   prox!(s1, ψ, mν∇fk, ν₁)
@@ -372,7 +381,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.σ = σk
+    solver.subpb.model.data.σ = σk
     isa(solver.subsolver, R2DHSolver) && (solver.subsolver.D.d[1] = 1/ν₁)
     if isa(solver.subsolver, R2Solver) #FIXME
       solve!(
@@ -458,9 +467,9 @@ function SolverCore.solve!(
       end
 
       if opnorm_maxiter ≤ 0
-        λmax, found_λ = opnorm(solver.subpb.model.B)
+        λmax, found_λ = opnorm(solver.subpb.model.data.H)
       else
-        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+        λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
       end
       found_λ || error("operator norm computation failed")
       set_step_status!(stats, :accepted)
diff --git a/src/R2NModel.jl b/src/R2NModel.jl
deleted file mode 100644
index 3ce4d779..00000000
--- a/src/R2NModel.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-export R2NModel
-
-@doc raw"""
-    R2NModel(B, ∇f, v, σ, x0)
-
-Given the unconstrained optimization problem:
-```math
-\min f(x),
-```
-this model represents the smooth R2N subproblem:
-```math
-\min_s \ ∇f^T s + \tfrac{1}{2} s^T B s + \tfrac{1}{2} σ \|s\|^2
-```
-where `B` is either an approximation of the Hessian of `f` or the Hessian itself and `∇f` represents the gradient of `f` at `x0`.
-`σ > 0` is a regularization parameter and `v` is a vector of the same size as `x0` used for intermediary computations.
-"""
-mutable struct R2NModel{T <: Real, V <: AbstractVector{T}, G <: AbstractLinearOperator{T}} <:
-               AbstractNLPModel{T, V}
-  B::G
-  ∇f::V
-  v::V
-  σ::T
-  meta::NLPModelMeta{T, V}
-  counters::Counters
-end
-
-function R2NModel(B::G, ∇f::V, σ::T, x0::V) where {T, V, G}
-  @assert length(x0) == length(∇f)
-  meta = NLPModelMeta(
-    length(∇f),
-    x0 = x0, # Perhaps we should add lvar and uvar as well here.
-  )
-  v = similar(x0)
-  return R2NModel(B::G, ∇f::V, v::V, σ::T, meta, Counters())
-end
-
-function NLPModels.obj(nlp::R2NModel, x::AbstractVector; skip_sigma = false)
-  @lencheck nlp.meta.nvar x
-  increment!(nlp, :neval_obj)
-  mul!(nlp.v, nlp.B, x)
-  skip_sigma && return dot(nlp.v, x)/2 + dot(nlp.∇f, x)
-  return dot(nlp.v, x)/2 + dot(nlp.∇f, x) + nlp.σ * dot(x, x) / 2
-end
-
-function NLPModels.grad!(nlp::R2NModel, x::AbstractVector, g::AbstractVector)
-  @lencheck nlp.meta.nvar x
-  @lencheck nlp.meta.nvar g
-  increment!(nlp, :neval_grad)
-  mul!(g, nlp.B, x)
-  g .+= nlp.∇f
-  g .+= nlp.σ .* x
-  return g
-end
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index e65d167c..de226752 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -11,6 +11,7 @@ using LinearOperators,
   ManualNLPModels,
   NLPModels,
   NLPModelsModifiers,
+  QuadraticModels,
   RegularizedProblems,
   ShiftedProximalOperators,
   SolverCore
@@ -45,7 +46,6 @@ include("LMModel.jl")
 include("LM_alg.jl")
 include("LMTR_alg.jl")
 include("R2DH.jl")
-include("R2NModel.jl")
 include("R2N.jl")
 include("AL_alg.jl")
 
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 7d3d9bbe..7ebec003 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -19,6 +19,7 @@ mutable struct TRSolver{
   xkn::V
   s::V
   v0::V
+  v1::V
   has_bnds::Bool
   l_bound::V
   u_bound::V
@@ -61,6 +62,7 @@ function TRSolver(
 
   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 ?
@@ -68,7 +70,7 @@ function TRSolver(
     shifted(reg_nlp.h, xk, T(1), χ)
 
   Bk = hess_op(reg_nlp, xk)
-  sub_nlp = R2NModel(Bk, ∇fk, zero(T), x0) #FIXME 
+  sub_nlp = QuadraticModel(∇fk, Bk, x0 = x0)
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
   subsolver = subsolver(subpb)
@@ -83,6 +85,7 @@ function TRSolver(
     xkn,
     s,
     v0,
+    v1,
     has_bnds,
     l_bound,
     u_bound,
@@ -307,9 +310,9 @@ function SolverCore.solve!(
   found_λ = true
 
   if opnorm_maxiter ≤ 0
-    λmax, found_λ = opnorm(solver.subpb.model.B)
+    λmax, found_λ = opnorm(solver.subpb.model.data.H)
   else
-    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+    λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
   end
   found_λ || error("operator norm computation failed")
 
@@ -468,9 +471,9 @@ function SolverCore.solve!(
       end
 
       if opnorm_maxiter ≤ 0
-        λmax, found_λ = opnorm(solver.subpb.model.B)
+        λmax, found_λ = opnorm(solver.subpb.model.data.H)
       else
-        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+        λmax = power_method!(solver.subpb.model.data.H, solver.v0, solver.v1, opnorm_maxiter)
       end
       found_λ || error("operator norm computation failed")