From fb2337660f1264c2b22391b7d54c1486cec87575 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Wed, 30 Oct 2024 23:52:27 -0400 Subject: [PATCH 1/8] Reuse Arpack in order to compute norm of Bk --- Project.toml | 4 +-- src/LMTR_alg.jl | 6 +++-- src/LM_alg.jl | 6 +++-- src/RegularizedOptimization.jl | 2 +- src/TR_alg.jl | 6 +++-- src/utils.jl | 47 ++++++++++++++++++++++++++++++++-- 6 files changed, 60 insertions(+), 11 deletions(-) diff --git a/Project.toml b/Project.toml index 97784bb0..98f6f703 100644 --- a/Project.toml +++ b/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] diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl index c4cfa838..6eff53c1 100644 --- a/src/LMTR_alg.jl +++ b/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 diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 56bc5356..59a3e579 100644 --- a/src/LM_alg.jl +++ b/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 diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl index f6953857..aab4a8a2 100644 --- a/src/RegularizedOptimization.jl +++ b/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, diff --git a/src/TR_alg.jl b/src/TR_alg.jl index 700b22c1..cd12aa0c 100644 --- a/src/TR_alg.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index f0853b7d..bc17dfa3 100644 --- a/src/utils.jl +++ b/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!( From 51853a385601507518cd5dc27c58931ceb37fbd6 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH <81633807+MohamedLaghdafHABIBOULLAH@users.noreply.github.com> Date: Thu, 31 Oct 2024 15:17:15 -0400 Subject: [PATCH 2/8] Update src/utils.jl Co-authored-by: Dominique --- src/utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils.jl b/src/utils.jl index bc17dfa3..61634abf 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -29,6 +29,7 @@ function opnorm_eig(B; max_attempts::Int = 3) end return λ, success end + function opnorm_svd(J; max_attempts::Int = 3) have_svd = false attempt = 0 From a35feb953be79367e0aefc0530137cda252d4e03 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Thu, 31 Oct 2024 15:26:08 -0400 Subject: [PATCH 3/8] Update RegularizedOptimization.jl and utils.jl Use alphabetical order for external dependencies and remove unused variables. --- src/RegularizedOptimization.jl | 2 +- src/utils.jl | 8 ++------ 2 files changed, 3 insertions(+), 7 deletions(-) diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl index aab4a8a2..53365e55 100644 --- a/src/RegularizedOptimization.jl +++ b/src/RegularizedOptimization.jl @@ -4,7 +4,7 @@ module RegularizedOptimization using LinearAlgebra, Logging, Printf # external dependencies -using ProximalOperators, Arpack +using Arpack, ProximalOperators # dependencies from us using LinearOperators, diff --git a/src/utils.jl b/src/utils.jl index 61634abf..be0f6d87 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -12,7 +12,6 @@ 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) @@ -22,19 +21,17 @@ function opnorm_eig(B; max_attempts::Int = 3) have_eig = nconv == 1 if (have_eig) λ = abs(d[1]) - success = true else ncv = min(2 * ncv, n) end end - return λ, success + return λ, have_eig 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 @@ -44,12 +41,11 @@ function opnorm_svd(J; max_attempts::Int = 3) have_svd = nconv == 1 if (have_svd) σ = maximum(s.S) - success = true else ncv = min(2 * ncv, n) end end - return σ, success + return σ, have_svd end ShiftedProximalOperators.iprox!( From e49f82b4e56552e5f9d03f40486dbfd9d264cd3f Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Mon, 4 Nov 2024 13:23:14 -0500 Subject: [PATCH 4/8] Update utils.jl 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. --- src/utils.jl | 70 ++++++++++++++++++++++++++++++++++++++-------------- 1 file changed, 51 insertions(+), 19 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index be0f6d87..bdc06359 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -15,16 +15,32 @@ function opnorm_eig(B; max_attempts::Int = 3) 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]) - else - ncv = min(2 * ncv, n) - end + + while !(have_eig || attempt >= max_attempts) + attempt += 1 + try + # Perform eigendecomposition + 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)) + @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 @@ -32,19 +48,35 @@ function opnorm_svd(J; max_attempts::Int = 3) have_svd = false attempt = 0 σ = zero(eltype(J)) - n = min(size(J)...) + n = min(size(J)...) # Minimum dimension of the matrix 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) - else - ncv = min(2 * ncv, n) - end + + while !(have_svd || attempt >= max_attempts) + attempt += 1 + try + # Perform singular value decomposition + 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)) + @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 From f36816474295017756845a7c77c461afc7dd3683 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH <81633807+MohamedLaghdafHABIBOULLAH@users.noreply.github.com> Date: Thu, 20 Feb 2025 18:56:31 -0500 Subject: [PATCH 5/8] Apply suggestions from code review Co-authored-by: Dominique --- src/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index bdc06359..6969ce52 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -19,7 +19,7 @@ function opnorm_eig(B; max_attempts::Int = 3) while !(have_eig || attempt >= max_attempts) attempt += 1 try - # Perform eigendecomposition + # 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 @@ -55,7 +55,7 @@ function opnorm_svd(J; max_attempts::Int = 3) while !(have_svd || attempt >= max_attempts) attempt += 1 try - # Perform singular value decomposition + # 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 From d763256e02a579f9e90eb3ff66eee1737a18e96a Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Thu, 20 Feb 2025 19:02:15 -0500 Subject: [PATCH 6/8] rethrow error if ncv = n --- src/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 6969ce52..7c5f351e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -32,7 +32,7 @@ function opnorm_eig(B; max_attempts::Int = 3) ncv = min(2 * ncv, n) end catch e - if occursin("XYAUPD_Exception", string(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 @@ -68,7 +68,7 @@ function opnorm_svd(J; max_attempts::Int = 3) ncv = min(2 * ncv, n) end catch e - if occursin("XYAUPD_Exception", string(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 From 3269c2804ba3c92e1d066ba7971fa8bd242edd41 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Wed, 26 Feb 2025 02:03:39 -0500 Subject: [PATCH 7/8] compute norm(Bk) within R2N --- src/R2N.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/R2N.jl b/src/R2N.jl index ba27f963..aba356a1 100644 --- a/src/R2N.jl +++ b/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 From 03591a4de1fd2b24c24035c6a42dae3ea52c63bb Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH <81633807+MohamedLaghdafHABIBOULLAH@users.noreply.github.com> Date: Wed, 26 Feb 2025 02:03:52 -0500 Subject: [PATCH 8/8] Update src/utils.jl Co-authored-by: Dominique --- src/utils.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/utils.jl b/src/utils.jl index 7c5f351e..a02e6043 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -8,6 +8,7 @@ function LinearAlgebra.opnorm(B; kwargs...) 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