Reuse Arpack in order to compute norm of Bk (#159)

Enhance robustness of Arpack usage by catching NCV-related errors

Catch the XYAUPD_Exception error in Arpack calls and dynamically increase NCV as needed to handle convergence issues. This modification addresses instability in LSR1 computations due to insufficient NCV.

Co-authored-by: Dominique <[email protected]>

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/R2N.jl

@@ -132,7 +132,8 @@ function R2N(
   quasiNewtTest = isa(f, QuasiNewtonModel)
   Bk = hess_op(f, xk)
 
-  λmax = opnorm(Bk)
+  λmax, found_λ = opnorm(Bk)
+  found_λ || error("operator norm computation failed")
   νInv = (1 + θ) * (σk + λmax)
   sqrt_ξ1_νInv = one(R)
 
@@ -247,7 +248,8 @@ function R2N(
         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/RegularizedOptimization.jl

@@ -4,7 +4,7 @@ module RegularizedOptimization
 using LinearAlgebra, Logging, Printf
 
 # external dependencies
-using ProximalOperators, TSVD
+using Arpack, ProximalOperators
 
 # 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,81 @@ 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))
+  n = size(B, 1)
+  nev = 1
+  ncv = max(20, 2 * nev + 1)
+
+  while !(have_eig || attempt >= max_attempts)
+      attempt += 1
+      try
+          # Estimate largest eigenvalue in absolute value
+          d, nconv, niter, nmult, resid = eigs(B; nev = nev, ncv = ncv, which = :LM, ritzvec = false, check = 1)
+          
+          # Check if eigenvalue has converged
+          have_eig = nconv == 1
+          if have_eig
+              λ = abs(d[1])  # Take absolute value of the largest eigenvalue
+              break  # Exit loop if successful
+          else
+              # Increase NCV for the next attempt if convergence wasn't achieved
+              ncv = min(2 * ncv, n)
+          end
+      catch e
+          if occursin("XYAUPD_Exception", string(e)) && ncv < n
+              @warn "Arpack error: $e. Increasing NCV to $ncv and retrying."
+              ncv = min(2 * ncv, n)  # Increase NCV but don't exceed matrix size
+          else
+              rethrow(e)  # Re-raise if it's a different error
+          end
+      end
+  end
+
+  return λ, have_eig
+end
+
+function opnorm_svd(J; max_attempts::Int = 3)
+  have_svd = false
+  attempt = 0
+  σ = zero(eltype(J))
+  n = min(size(J)...)  # Minimum dimension of the matrix
+  nsv = 1
+  ncv = 10
+
+  while !(have_svd || attempt >= max_attempts)
+      attempt += 1
+      try
+          # Estimate largest singular value
+          s, nconv, niter, nmult, resid = svds(J; nsv = nsv, ncv = ncv, ritzvec = false, check = 1)
+          
+          # Check if singular value has converged
+          have_svd = nconv >= 1
+          if have_svd
+              σ = maximum(s.S)  # Take the largest singular value
+              break  # Exit loop if successful
+          else
+              # Increase NCV for the next attempt if convergence wasn't achieved
+              ncv = min(2 * ncv, n)
+          end
+      catch e
+          if occursin("XYAUPD_Exception", string(e)) && ncv < n
+              @warn "Arpack error: $e. Increasing NCV to $ncv and retrying."
+              ncv = min(2 * ncv, n)  # Increase NCV but don't exceed matrix size
+          else
+              rethrow(e)  # Re-raise if it's a different error
+          end
+      end
+  end
+
+  return σ, have_svd
 end
 
 ShiftedProximalOperators.iprox!(