diff --git a/src/R2N.jl b/src/R2N.jl
index acb05d94..2250547f 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -35,6 +35,7 @@ function R2NSolver(
   reg_nlp::AbstractRegularizedNLPModel{T, V};
   subsolver = R2Solver,
   m_monotone::Int = 1,
+  sparse = false
 ) where {T, V}
   x0 = reg_nlp.model.meta.x0
   l_bound = reg_nlp.model.meta.lvar
@@ -69,7 +70,8 @@ function R2NSolver(
     has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) :
     shifted(reg_nlp.h, xk)
 
-  Bk = hess_op(reg_nlp, xk)
+  sparse = isa(reg_nlp.model, QuasiNewtonModel) ? false : sparse
+  Bk = sparse ? hess(reg_nlp, xk) : hess_op(reg_nlp, xk)
   sub_nlp = QuadraticModel(∇fk, Bk, σ = T(1), x0 = x0)
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
@@ -464,6 +466,8 @@ function SolverCore.solve!(
         qn_update_y!(nlp, solver, stats)
         push!(nlp, s, solver.y)
         qn_copy!(nlp, solver, stats)
+      elseif isa(solver.subpb.model.data.H, AbstractMatrix)
+        hess_coord!(reg_nlp, xk, solver.subpb.model.data.H.nzval)
       end
 
       if opnorm_maxiter ≤ 0
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 7ebec003..20347f3b 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -36,6 +36,7 @@ function TRSolver(
   χ::X = NormLinf(one(T)),
   subsolver = R2Solver,
   m_monotone::Int = 1,
+  sparse::Bool = false,
 ) where {T, V, X}
   x0 = reg_nlp.model.meta.x0
   l_bound = reg_nlp.model.meta.lvar
@@ -69,7 +70,8 @@ function TRSolver(
     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)
+  sparse = isa(reg_nlp.model, QuasiNewtonModel) ? false : sparse
+  Bk = sparse ? hess(reg_nlp, xk) : hess_op(reg_nlp, xk)
   sub_nlp = QuadraticModel(∇fk, Bk, x0 = x0)
   subpb = RegularizedNLPModel(sub_nlp, ψ)
   substats = RegularizedExecutionStats(subpb)
@@ -468,6 +470,8 @@ function SolverCore.solve!(
       if quasiNewtTest
         @. ∇fk⁻ = ∇fk - ∇fk⁻
         push!(nlp, s, ∇fk⁻) # update QN operator
+      elseif isa(solver.subpb.model.data.H, AbstractMatrix)
+        hess_coord!(reg_nlp, xk, solver.subpb.model.data.H.nzval)
       end
 
       if opnorm_maxiter ≤ 0