diff --git a/Project.toml b/Project.toml
index 8208c222..02f46a80 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ author = ["Robert Baraldi <rbaraldi@uw.edu> and Dominique Orban <dominique.orban
 version = "0.1.0"
 
 [deps]
-Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
@@ -13,16 +12,17 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
+TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 
 [compat]
-Arpack = "0.5"
 NLPModels = "0.19"
 NLPModelsModifiers = "0.6"
 ProximalOperators = "0.15"
 RegularizedProblems = "0.1"
 ShiftedProximalOperators = "0.1"
 SolverCore = "0.3.0"
-julia = "^1.3.0"
+TSVD = "0.4"
+julia = "^1.6.0"
 
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl
index dd44136d..ee4609f7 100644
--- a/src/LMTR_alg.jl
+++ b/src/LMTR_alg.jl
@@ -110,8 +110,7 @@ function LMTR(
   JdFk = similar(Fk)   # temporary storage
   Jt_Fk = similar(∇fk)   # temporary storage
 
-  σmax, found_σ = opnorm(Jk)
-  found_σ || error("operator norm computation failed")
+  σmax = opnorm(Jk)
   νInv = (1 + θ) * σmax^2  # ‖J'J‖ = ‖J‖²
 
   mν∇fk = -∇fk / νInv
@@ -217,8 +216,7 @@ function LMTR(
       shift!(ψ, xk)
       Jk = jac_op_residual(nls, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
-      σmax, found_σ = opnorm(Jk)
-      found_σ || error("operator norm computation failed")
+      σmax = opnorm(Jk)
       νInv = (1 + θ) * σmax^2  # ‖J'J‖ = ‖J‖²
       @. mν∇fk = -∇fk / νInv
     end
diff --git a/src/LM_alg.jl b/src/LM_alg.jl
index 24b4612a..e9e6b5ed 100644
--- a/src/LM_alg.jl
+++ b/src/LM_alg.jl
@@ -106,8 +106,7 @@ function LM(
   JdFk = similar(Fk)   # temporary storage
   Jt_Fk = similar(∇fk)
 
-  σmax, found_σ = opnorm(Jk)
-  found_σ || error("operator norm computation failed")
+  σmax = opnorm(Jk)
   νInv = (1 + θ) * (σmax^2 + σk)  # ‖J'J + σₖ I‖ = ‖J‖² + σₖ
 
   s = zero(xk)
@@ -217,8 +216,7 @@ function LM(
       Jk = jac_op_residual(nls, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
 
-      σmax, found_σ = opnorm(Jk)
-      found_σ || error("operator norm computation failed")
+      σmax = opnorm(Jk)
       νInv = (1 + θ) * (σmax^2 + σk)  # ‖J'J + σₖ I‖ = ‖J‖² + σₖ
 
       Complex_hist[k] += 1
diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl
index 4854dc91..d71d32d2 100644
--- a/src/RegularizedOptimization.jl
+++ b/src/RegularizedOptimization.jl
@@ -4,7 +4,7 @@ module RegularizedOptimization
 using LinearAlgebra, Logging, Printf
 
 # external dependencies
-using Arpack, ProximalOperators
+using ProximalOperators, TSVD
 
 # dependencies from us
 using NLPModels, NLPModelsModifiers, ShiftedProximalOperators, SolverCore
diff --git a/src/TR_alg.jl b/src/TR_alg.jl
index c1195a32..829983c9 100644
--- a/src/TR_alg.jl
+++ b/src/TR_alg.jl
@@ -114,8 +114,7 @@ function TR(
   quasiNewtTest = isa(f, QuasiNewtonModel)
   Bk = hess_op(f, xk)
 
-  λmax, found_λ = opnorm(Bk)
-  found_λ || error("operator norm computation failed")
+  λmax = opnorm(Bk)
   νInv = (1 + θ) * λmax
 
   optimal = false
@@ -207,8 +206,7 @@ function TR(
         push!(f, s, ∇fk - ∇fk⁻)
       end
       Bk = hess_op(f, xk)
-      λmax, found_λ = opnorm(Bk)
-      found_λ || error("operator norm computation failed")
+      λmax = opnorm(Bk)
       νInv = (1 + θ) * λmax
       ∇fk⁻ .= ∇fk
     end
diff --git a/src/utils.jl b/src/utils.jl
index 42690811..f2913a75 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,51 +1,5 @@
 # use Arpack to obtain largest eigenvalue in magnitude with a minimum of robustness
 function LinearAlgebra.opnorm(B; kwargs...)
-  m, n = size(B)
-  opnorm_fcn = m == n ? opnorm_eig : opnorm_svd
-  return opnorm_fcn(B; kwargs...)
-end
-
-function opnorm_eig(B; max_attempts::Int = 3)
-  have_eig = false
-  attempt = 0
-  λ = zero(eltype(B))
-  success = false
-  n = size(B, 1)
-  nev = 1
-  ncv = max(20, 2 * nev + 1)
-  while !(have_eig || attempt > max_attempts)
-    attempt += 1
-    (d, nconv, niter, nmult, resid) =
-      eigs(B; nev = nev, ncv = ncv, which = :LM, ritzvec = false, check = 1)
-    have_eig = nconv == 1
-    if (have_eig)
-      λ = abs(d[1])
-      success = true
-    else
-      ncv = min(2 * ncv, n)
-    end
-  end
-  return λ, success
-end
-
-function opnorm_svd(J; max_attempts::Int = 3)
-  have_svd = false
-  attempt = 0
-  σ = zero(eltype(J))
-  success = false
-  n = min(size(J)...)
-  nsv = 1
-  ncv = 10
-  while !(have_svd || attempt > max_attempts)
-    attempt += 1
-    (s, nconv, niter, nmult, resid) = svds(J, nsv = nsv, ncv = ncv, ritzvec = false, check = 1)
-    have_svd = nconv == 1
-    if (have_svd)
-      σ = maximum(s.S)
-      success = true
-    else
-      ncv = min(2 * ncv, n)
-    end
-  end
-  return σ, success
+  _, s, _ = tsvd(B)
+  return s[1]
 end