diff --git a/src/R2N.jl b/src/R2N.jl
index cb6c066f..a6c52f5e 100644
--- a/src/R2N.jl
+++ b/src/R2N.jl
@@ -18,6 +18,7 @@ mutable struct R2NSolver{
   xkn::V
   s::V
   s1::V
+  v0::V
   has_bnds::Bool
   l_bound::V
   u_bound::V
@@ -46,6 +47,10 @@ function R2NSolver(
   xkn = similar(x0)
   s = similar(x0)
   s1 = similar(x0)
+
+  v0 = [(-1.0)^i for i in 0:(reg_nlp.model.meta.nvar-1)]
+  v0 ./= sqrt(reg_nlp.model.meta.nvar)
+
   has_bnds = any(l_bound .!= T(-Inf)) || any(u_bound .!= T(Inf))
   if has_bnds
     l_bound_m_x = similar(xk)
@@ -78,6 +83,7 @@ function R2NSolver(
     xkn,
     s,
     s1,
+    v0,
     has_bnds,
     l_bound,
     u_bound,
@@ -133,6 +139,7 @@ For advanced usage, first define a solver "R2NSolver" to preallocate the memory
 - `η2::T = T(0.9)`: very successful iteration threshold;
 - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful;
 - `θ::T = 1/(1 + eps(T)^(1 / 5))`: is the model decrease fraction with respect to the decrease of the Cauchy model;
+- `opnorm_maxiter::Int = 5`: how many iterations of the power method to use to compute the operator norm of Bₖ. If a negative number is provided, then Arpack is used instead;
 - `m_monotone::Int = 1`: monotonicity parameter. By default, R2N is monotone but the non-monotone variant will be used if `m_monotone > 1`;
 - `sub_kwargs::NamedTuple = NamedTuple()`: a named tuple containing the keyword arguments to be sent to the subsolver. The solver will fail if invalid keyword arguments are provided to the subsolver. For example, if the subsolver is `R2Solver`, you can pass `sub_kwargs = (max_iter = 100, σmin = 1e-6,)`.
 
@@ -160,7 +167,6 @@ function R2N(
   selected = pop!(kwargs_dict, :selected, 1:(nlp.meta.nvar))
   x0 = pop!(kwargs_dict, :x0, nlp.meta.x0)
   reg_nlp = RegularizedNLPModel(nlp, h, selected)
-  sub_kwargs = pop!(kwargs_dict, :sub_kwargs, NamedTuple())
   return R2N(
     reg_nlp,
     x = x0,
@@ -174,8 +180,7 @@ function R2N(
     σk = options.σk,
     η1 = options.η1,
     η2 = options.η2,
-    γ = options.γ,
-    sub_kwargs = sub_kwargs;
+    γ = options.γ;
     kwargs_dict...,
   )
 end
@@ -212,6 +217,7 @@ function SolverCore.solve!(
   γ::T = T(3),
   β::T = 1 / eps(T),
   θ::T = 1/(1 + eps(T)^(1 / 5)),
+  opnorm_maxiter::Int = 5,
   sub_kwargs::NamedTuple = NamedTuple(),
 ) where {T, V, G}
   reset!(stats)
@@ -285,9 +291,14 @@ function SolverCore.solve!(
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
   λmax::T = T(1)
+  found_λ = true
   solver.subpb.model.B = hess_op(nlp, xk)
 
-  λmax, found_λ = opnorm(solver.subpb.model.B)
+  if opnorm_maxiter ≤ 0
+    λmax, found_λ = opnorm(solver.subpb.model.B)
+  else
+    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+  end
   found_λ || error("operator norm computation failed")
 
   ν₁ = θ / (λmax + σk)
@@ -437,7 +448,11 @@ function SolverCore.solve!(
       end
       solver.subpb.model.B = hess_op(nlp, xk)
 
-      λmax, found_λ = opnorm(solver.subpb.model.B)
+      if opnorm_maxiter ≤ 0
+        λmax, found_λ = opnorm(solver.subpb.model.B)
+      else
+        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+      end
       found_λ || error("operator norm computation failed")
     end
 
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index 6e8bb42b..aa9cca0d 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -18,6 +18,7 @@ mutable struct TRSolver{
   χ::N
   xkn::V
   s::V
+  v0::V
   has_bnds::Bool
   l_bound::V
   u_bound::V
@@ -58,6 +59,9 @@ function TRSolver(
 
   m_fh_hist = fill(T(-Inf), m_monotone - 1)
 
+  v0 = [(-1.0)^i for i in 0:(reg_nlp.model.meta.nvar-1)]
+  v0 ./= sqrt(reg_nlp.model.meta.nvar)
+
   ψ =
     has_bnds || subsolver == TRDHSolver ?
     shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) :
@@ -80,6 +84,7 @@ function TRSolver(
     χ,
     xkn,
     s,
+    v0,
     has_bnds,
     l_bound,
     u_bound,
@@ -135,6 +140,7 @@ For advanced usage, first define a solver "TRSolver" to preallocate the memory u
 - `η2::T = T(0.9)`: very successful iteration threshold;
 - `γ::T = T(3)`: trust-region radius parameter multiplier. Must satisfy `γ > 1`. The trust-region radius is updated as Δ := Δ*γ when the iteration is very successful and Δ := Δ/γ when the iteration is unsuccessful;
 - `m_monotone::Int = 1`: monotonicity parameter. By default, TR is monotone but the non-monotone variant will be used if `m_monotone > 1`;
+- `opnorm_maxiter::Int = 5`: how many iterations of the power method to use to compute the operator norm of Bₖ. If a negative number is provided, then Arpack is used instead;
 - `χ::F =  NormLinf(1)`: norm used to define the trust-region;`
 - `subsolver::S = R2Solver`: subsolver used to solve the subproblem that appears at each iteration.
 - `sub_kwargs::NamedTuple = NamedTuple()`: a named tuple containing the keyword arguments to be sent to the subsolver. The solver will fail if invalid keyword arguments are provided to the subsolver. For example, if the subsolver is `R2Solver`, you can pass `sub_kwargs = (max_iter = 100, σmin = 1e-6,)`.
@@ -208,6 +214,7 @@ function SolverCore.solve!(
   η2::T = T(0.9),
   γ::T = T(3),
   sub_kwargs::NamedTuple = NamedTuple(),
+  opnorm_maxiter::Int = 5,
 ) where {T, G, V}
   reset!(stats)
 
@@ -286,9 +293,14 @@ function SolverCore.solve!(
   ∇fk⁻ .= ∇fk
 
   quasiNewtTest = isa(nlp, QuasiNewtonModel)
-  λmax = T(1)
+  λmax::T = T(1)
+  found_λ = true
 
-  λmax, found_λ = opnorm(solver.subpb.model.B)
+  if opnorm_maxiter ≤ 0
+    λmax, found_λ = opnorm(solver.subpb.model.B)
+  else
+    λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+  end
   found_λ || error("operator norm computation failed")
 
   ν₁ = α * Δk / (1 + λmax * (α * Δk + 1))
@@ -345,7 +357,6 @@ function SolverCore.solve!(
   callback(nlp, solver, stats)
 
   done = stats.status != :unknown
-
   while !done
     sub_atol = stats.iter == 0 ? 1e-5 : max(sub_atol, min(1e-2, sqrt_ξ1_νInv))
     ∆_effective = min(β * χ(s), Δk)
@@ -446,7 +457,11 @@ function SolverCore.solve!(
         push!(nlp, s, ∇fk⁻) # update QN operator
       end
 
-      λmax, found_λ = opnorm(solver.subpb.model.B)
+      if opnorm_maxiter ≤ 0
+        λmax, found_λ = opnorm(solver.subpb.model.B)
+      else
+        λmax = power_method!(solver.subpb.model.B, solver.v0, solver.subpb.model.v, opnorm_maxiter)
+      end
       found_λ || error("operator norm computation failed")
 
       ∇fk⁻ .= ∇fk
diff --git a/src/utils.jl b/src/utils.jl
index acc35592..a018abeb 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -2,6 +2,19 @@ export RegularizedExecutionStats
 
 import SolverCore.GenericExecutionStats
 
+function power_method!(B::M, v₀::S, v₁::S, max_iter::Int = 1) where{M, S}
+  @assert max_iter >= 1 "max_iter must be at least 1."
+  mul!(v₁, B, v₀)
+  normalize!(v₁) # v1 = B*v0 / ‖B*v0‖
+  for i = 2:max_iter
+    v₀ .= v₁ # v0 = v1
+    mul!(v₁, B, v₀)
+    normalize!(v₁)
+  end
+  mul!(v₁, B, v₀)
+  return abs(dot(v₀, v₁))
+end
+
 # use Arpack to obtain largest eigenvalue in magnitude with a minimum of robustness
 function LinearAlgebra.opnorm(B; kwargs...)
   m, n = size(B)
diff --git a/test/runtests.jl b/test/runtests.jl
index 4daaa1f5..071616c0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,6 +10,7 @@ using ADNLPModels,
   RegularizedOptimization,
   SolverCore
 
+Random.seed!(0)
 const global compound = 1
 const global nz = 10 * compound
 const global options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
diff --git a/test/test_allocs.jl b/test/test_allocs.jl
index 795d4528..0bc22dee 100644
--- a/test/test_allocs.jl
+++ b/test/test_allocs.jl
@@ -47,20 +47,18 @@ end
       (:R2DHSolver, "R2DH"),
       (:R2NSolver, "R2N"),
       (:TRDHSolver, "TRDH"),
-      (:TRSolver, "TR"),
     )
       @testset "$(solver_name)" begin
-        (solver_name == "R2N" || solver_name == "TR") && continue #FIXME
         reg_nlp = RegularizedNLPModel(LBFGSModel(bpdn), h)
         solver = eval(solver)(reg_nlp)
         stats = RegularizedExecutionStats(reg_nlp)
         solver_name == "R2" &&
           @test @wrappedallocs(solve!(solver, reg_nlp, stats, ν = 1.0, atol = 1e-6, rtol = 1e-6)) ==
                 0
-        solver_name == "R2DH" && @test @wrappedallocs(
+        (solver_name == "R2DH" || solver_name == "R2N") && @test @wrappedallocs(
           solve!(solver, reg_nlp, stats, σk = 1.0, atol = 1e-6, rtol = 1e-6)
         ) == 0
-        solver_name == "TRDH" &&
+        (solver_name == "TRDH") &&
           @test @wrappedallocs(solve!(solver, reg_nlp, stats, atol = 1e-6, rtol = 1e-6)) == 0
         @test stats.status == :first_order
       end