Reuse Arpack in order to compute norm of Bk

Author: MohamedLaghdafHABIBOULLAH

Project.toml

@@ -15,7 +15,7 @@ ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
 RegularizedProblems = "ea076b23-609f-44d2-bb12-a4ae45328278"
 ShiftedProximalOperators = "d4fd37fa-580c-4e43-9b30-361c21aae263"
 SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
-TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
+Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
 
 [compat]
 LinearOperators = "2.9"
@@ -26,7 +26,7 @@ ProximalOperators = "0.15"
 RegularizedProblems = "0.1.1"
 ShiftedProximalOperators = "0.2"
 SolverCore = "0.3.0"
-TSVD = "0.4"
+Arpack = "0.5"
 julia = "^1.6.0"
 
 [extras]

src/LMTR_alg.jl

@@ -130,7 +130,8 @@ function LMTR(
   JdFk = similar(Fk)   # temporary storage
   Jt_Fk = similar(∇fk)   # temporary storage
 
-  σmax = opnorm(Jk)
+  σmax, found_σ = opnorm(Jk)
+  found_σ || error("operator norm computation failed")
   νInv = (1 + θ) * σmax^2  # ‖J'J‖ = ‖J‖²
 
   mν∇fk = -∇fk / νInv
@@ -252,7 +253,8 @@ function LMTR(
       shift!(ψ, xk)
       Jk = jac_op_residual(nls, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
-      σmax = opnorm(Jk)
+      σmax, found_σ = opnorm(Jk)
+      found_σ || error("operator norm computation failed")
       νInv = (1 + θ) * σmax^2  # ‖J'J‖ = ‖J‖²
       @. mν∇fk = -∇fk / νInv
     end

src/LM_alg.jl

@@ -123,7 +123,8 @@ function LM(
   JdFk = similar(Fk)   # temporary storage
   Jt_Fk = similar(∇fk)
 
-  σmax = opnorm(Jk)
+  σmax, found_σ = opnorm(Jk)
+  found_σ || error("operator norm computation failed")
   νInv = (1 + θ) * (σmax^2 + σk)  # ‖J'J + σₖ I‖ = ‖J‖² + σₖ
 
   s = zero(xk)
@@ -243,7 +244,8 @@ function LM(
       Jk = jac_op_residual(nls, xk)
       jtprod_residual!(nls, xk, Fk, ∇fk)
 
-      σmax = opnorm(Jk)
+      σmax, found_σ = opnorm(Jk)
+      found_σ || error("operator norm computation failed")
 
       Complex_hist[k] += 1
     end

src/RegularizedOptimization.jl

@@ -4,7 +4,7 @@ module RegularizedOptimization
 using LinearAlgebra, Logging, Printf
 
 # external dependencies
-using ProximalOperators, TSVD
+using ProximalOperators, Arpack
 
 # dependencies from us
 using LinearOperators,

src/TR_alg.jl

@@ -131,7 +131,8 @@ function TR(
   quasiNewtTest = isa(f, QuasiNewtonModel)
   Bk = hess_op(f, xk)
 
-  λmax = opnorm(Bk)
+  λmax, found_λ = opnorm(Bk)
+  found_λ || error("operator norm computation failed")
   α⁻¹Δ⁻¹ = 1 / (α * Δk)
   ν = 1 / (α⁻¹Δ⁻¹ + λmax * (α⁻¹Δ⁻¹ + 1))
   sqrt_ξ1_νInv = one(R)
@@ -241,7 +242,8 @@ function TR(
         push!(f, s, ∇fk - ∇fk⁻)
       end
       Bk = hess_op(f, xk)
-      λmax = opnorm(Bk)
+      λmax, found_λ = opnorm(Bk)
+      found_λ || error("operator norm computation failed")
       ∇fk⁻ .= ∇fk
     end
 

src/utils.jl

@@ -4,8 +4,51 @@ import SolverCore.GenericExecutionStats
 
 # use Arpack to obtain largest eigenvalue in magnitude with a minimum of robustness
 function LinearAlgebra.opnorm(B; kwargs...)
-  _, s, _ = tsvd(B)
-  return s[1]
+  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
 end
 
 ShiftedProximalOperators.iprox!(