From 2130fc9ab7f9737cf1c8ae79f2e7b77df7b3d5b5 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Fri, 27 Sep 2024 17:44:37 -0400 Subject: [PATCH 1/6] update LM solver with possibility of using R2DH as subsolver and with different subsolver criterion stop --- src/LM_alg.jl | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 59a3e579..2292b930 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -47,6 +47,7 @@ function LM( subsolver = R2, subsolver_options = ROSolverOptions(ϵa = options.ϵa), selected::AbstractVector{<:Integer} = 1:(nls.meta.nvar), + nonlinear::Bool = true, ) where {H} start_time = time() elapsed_time = 0.0 @@ -62,6 +63,7 @@ function LM( γ = options.γ θ = options.θ σmin = options.σmin + σk = options.σk # store initial values of the subsolver_options fields that will be modified ν_subsolver = subsolver_options.ν @@ -85,7 +87,6 @@ function LM( end # initialize parameters - σk = max(1 / options.ν, σmin) xk = copy(x0) hk = h(xk[selected]) if hk == Inf @@ -104,13 +105,15 @@ function LM( k = 0 Fobj_hist = zeros(maxIter) Hobj_hist = zeros(maxIter) + R = eltype(xk) + sqrt_ξ1_νInv = one(R) Complex_hist = zeros(Int, maxIter) Grad_hist = zeros(Int, maxIter) Resid_hist = zeros(Int, maxIter) if verbose > 0 #! format: off - @info @sprintf "%6s %8s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "inner" "f(x)" "h(x)" "√ξ1" "√ξ" "ρ" "σ" "‖x‖" "‖s‖" "‖Jₖ‖²" "reg" + @info @sprintf "%6s %8s %8s %8s %7s %7s %8s %7s %7s %7s %7s %1s" "outer" "inner" "f(x)" "h(x)" "√(ξ1/ν)" "√ξ" "ρ" "σ" "‖x‖" "‖s‖" "‖Jₖ‖²" "reg" #! format: on end @@ -177,22 +180,24 @@ function LM( prox!(s, ψ, ∇fk, ν) ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps() # TODO: isn't mk(s) returned by subsolver? ξ1 > 0 || error("LM: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)") + sqrt_ξ1_νInv = ξ1 > 0 ? sqrt(ξ1 * νInv) : 10. if ξ1 ≥ 0 && k == 1 - ϵ_increment = ϵr * sqrt(ξ1) + ϵ_increment = ϵr * sqrt_ξ1_νInv ϵ += ϵ_increment # make stopping test absolute and relative ϵ_subsolver += ϵ_increment end - if sqrt(ξ1) < ϵ + if sqrt_ξ1_νInv < ϵ # the current xk is approximately first-order stationary optimal = true continue end - subsolver_options.ϵa = k == 1 ? 1.0e-1 : max(ϵ_subsolver, min(1.0e-2, ξ1 / 10)) + # subsolver_options.ϵa = k == 1 ? 1.0e-1 : max(ϵ_subsolver, min(1.0e-2, ξ1 / 10)) + subsolver_options.ϵa = k == 1 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5) , sqrt_ξ1_νInv * 1e-3) # 1.0e-5 default subsolver_options.ν = ν - subsolver_args = subsolver == TRDH ? (SpectralGradient(1 / ν, nls.meta.nvar),) : () + subsolver_args = subsolver == R2DH ? (SpectralGradient(νInv, nls.meta.nvar),) : () @debug "setting inner stopping tolerance to" subsolver_options.optTol s, iter, _ = with_logger(subsolver_logger) do subsolver(φ, ∇φ!, ψ, subsolver_args..., subsolver_options, s) @@ -222,7 +227,7 @@ function LM( if (verbose > 0) && (k % ptf == 0) #! format: off - @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt(ξ1) sqrt(ξ) ρk σk norm(xk) norm(s) νInv σ_stat + @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8.1e %7.1e %7.1e %7.1e %7.1e %1s" k iter fk hk sqrt_ξ1_νInv sqrt(ξ) ρk σk norm(xk) norm(s) νInv σ_stat #! format: off end @@ -244,8 +249,11 @@ function LM( Jk = jac_op_residual(nls, xk) jtprod_residual!(nls, xk, Fk, ∇fk) - σmax, found_σ = opnorm(Jk) - found_σ || error("operator norm computation failed") + # update opnorm if not linear least squares + if nonlinear == true + σmax, found_σ = opnorm(Jk) + found_σ || error("operator norm computation failed") + end Complex_hist[k] += 1 end @@ -254,6 +262,7 @@ function LM( σk = σk * γ end νInv = (1 + θ) * (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ + tired = k ≥ maxIter || elapsed_time > maxTime end @@ -262,9 +271,9 @@ function LM( @info @sprintf "%6d %8s %8.1e %8.1e" k "" fk hk elseif optimal #! format: off - @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k 1 fk hk sqrt(ξ1) sqrt(ξ1) "" σk norm(xk) norm(s) νInv + @info @sprintf "%6d %8d %8.1e %8.1e %7.1e %7.1e %8s %7.1e %7.1e %7.1e %7.1e" k 1 fk hk sqrt_ξ1_νInv sqrt(ξ1) "" σk norm(xk) norm(s) νInv #! format: on - @info "LM: terminating with √ξ1 = $(sqrt(ξ1))" + @info "LM: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" end end status = if optimal @@ -281,7 +290,7 @@ function LM( set_status!(stats, status) set_solution!(stats, xk) set_objective!(stats, fk + hk) - set_residuals!(stats, zero(eltype(xk)), ξ1 ≥ 0 ? sqrt(ξ1) : ξ1) + set_residuals!(stats, zero(eltype(xk)), ξ1 ≥ 0 ? sqrt_ξ1_νInv : ξ1) set_iter!(stats, k) set_time!(stats, elapsed_time) set_solver_specific!(stats, :Fhist, Fobj_hist[1:k]) From 825729d6c311d37807148df5e22ad54cfda369d4 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH <81633807+MohamedLaghdafHABIBOULLAH@users.noreply.github.com> Date: Thu, 20 Feb 2025 15:02:57 -0500 Subject: [PATCH 2/6] Apply suggestions from code review Co-authored-by: Dominique --- src/LM_alg.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 2292b930..49e39e60 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -262,7 +262,6 @@ function LM( σk = σk * γ end νInv = (1 + θ) * (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ - tired = k ≥ maxIter || elapsed_time > maxTime end From 435a0699631a4c0bfe06cc28896b31a4785a991c Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Thu, 20 Feb 2025 15:08:37 -0500 Subject: [PATCH 3/6] =?UTF-8?q?correct=20the=20definition=20of=20sqrt=5F?= =?UTF-8?q?=CE=BE1=5F=CE=BDInv?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/LM_alg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 49e39e60..a8da994c 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -180,7 +180,7 @@ function LM( prox!(s, ψ, ∇fk, ν) ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps() # TODO: isn't mk(s) returned by subsolver? ξ1 > 0 || error("LM: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - sqrt_ξ1_νInv = ξ1 > 0 ? sqrt(ξ1 * νInv) : 10. + sqrt_ξ1_νInv = sqrt(ξ1 * νInv) if ξ1 ≥ 0 && k == 1 ϵ_increment = ϵr * sqrt_ξ1_νInv From 1787811cb33fd2db9d42de57739a126a7cc58281 Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Thu, 20 Feb 2025 15:14:40 -0500 Subject: [PATCH 4/6] =?UTF-8?q?update=20=CE=BDInv=20according=20to=20Maxen?= =?UTF-8?q?ce=20Proposal?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/LM_alg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index a8da994c..9f963da2 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -261,7 +261,7 @@ function LM( if ρk < η1 || ρk == Inf σk = σk * γ end - νInv = (1 + θ) * (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ + νInv = (σmax^2 + σk) / θ # ‖J'J + σₖ I‖ = ‖J‖² + σₖ tired = k ≥ maxIter || elapsed_time > maxTime end From b1c8037428a64419e8c3612fc22d295da27557fc Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Sun, 23 Feb 2025 22:52:44 -0500 Subject: [PATCH 5/6] remove theta update, as it will be rebased later --- src/LM_alg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 9f963da2..a8da994c 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -261,7 +261,7 @@ function LM( if ρk < η1 || ρk == Inf σk = σk * γ end - νInv = (σmax^2 + σk) / θ # ‖J'J + σₖ I‖ = ‖J‖² + σₖ + νInv = (1 + θ) * (σmax^2 + σk) # ‖J'J + σₖ I‖ = ‖J‖² + σₖ tired = k ≥ maxIter || elapsed_time > maxTime end From 579cf06394ebe867c6bb8d45e672223475b3780b Mon Sep 17 00:00:00 2001 From: MohamedLaghdafHABIBOULLAH Date: Wed, 26 Feb 2025 04:47:15 -0500 Subject: [PATCH 6/6] declare sqrt_xi_1 as a local variable and add unit test for LM with R2DH as subsolver --- src/LM_alg.jl | 3 +- test/test_bounds.jl | 88 ++++++++++++++++++++++++++++++++------------- 2 files changed, 64 insertions(+), 27 deletions(-) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index a8da994c..163cc6fe 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -102,11 +102,10 @@ function LM( xkn = similar(xk) local ξ1 + local sqrt_ξ1_νInv k = 0 Fobj_hist = zeros(maxIter) Hobj_hist = zeros(maxIter) - R = eltype(xk) - sqrt_ξ1_νInv = one(R) Complex_hist = zeros(Int, maxIter) Grad_hist = zeros(Int, maxIter) Resid_hist = zeros(Int, maxIter) diff --git a/test/test_bounds.jl b/test/test_bounds.jl index d4ae95e8..2728820a 100644 --- a/test/test_bounds.jl +++ b/test/test_bounds.jl @@ -66,30 +66,68 @@ for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1")) @test h(out.solution) == out.solver_specific[:Hhist][end] @test out.status == :first_order end - @testset "bpdn-with-bounds-ls-$(solver_name)-$(h_name)-TRDH" begin - x0 = zeros(bpdn_nls2.meta.nvar) - args = solver_sym == :LM ? () : (NormLinf(1.0),) - @test has_bounds(bpdn_nls2) - out = solver( - bpdn_nls2, - h, - args..., - options, - x0 = x0, - subsolver = TRDH, - subsolver_options = subsolver_options, - ) - @test typeof(out.solution) == typeof(bpdn_nls2.meta.x0) - @test length(out.solution) == bpdn_nls2.meta.nvar - @test typeof(out.solver_specific[:Fhist]) == typeof(out.solution) - @test typeof(out.solver_specific[:Hhist]) == typeof(out.solution) - @test typeof(out.solver_specific[:SubsolverCounter]) == Array{Int, 1} - @test typeof(out.dual_feas) == eltype(out.solution) - @test length(out.solver_specific[:Fhist]) == length(out.solver_specific[:Hhist]) - @test length(out.solver_specific[:Fhist]) == length(out.solver_specific[:SubsolverCounter]) - @test obj(bpdn_nls2, out.solution) == out.solver_specific[:Fhist][end] - @test h(out.solution) == out.solver_specific[:Hhist][end] - @test out.status == :first_order - end + end +end + +# LMTR with TRDH as subsolver +for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1")) + @testset "bpdn-with-bounds-ls-LMTR-$(h_name)-TRDH" begin + x0 = zeros(bpdn_nls2.meta.nvar) + @test has_bounds(bpdn_nls2) + LMTR_out = LMTR( + bpdn_nls2, + h, + NormLinf(1.0), + options, + x0 = x0, + subsolver = TRDH, + subsolver_options = subsolver_options, + ) + @test typeof(LMTR_out.solution) == typeof(bpdn_nls2.meta.x0) + @test length(LMTR_out.solution) == bpdn_nls2.meta.nvar + @test typeof(LMTR_out.solver_specific[:Fhist]) == typeof(LMTR_out.solution) + @test typeof(LMTR_out.solver_specific[:Hhist]) == typeof(LMTR_out.solution) + @test typeof(LMTR_out.solver_specific[:SubsolverCounter]) == Array{Int, 1} + @test typeof(LMTR_out.dual_feas) == eltype(LMTR_out.solution) + @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:Hhist]) + @test length(LMTR_out.solver_specific[:Fhist]) == + length(LMTR_out.solver_specific[:SubsolverCounter]) + @test length(LMTR_out.solver_specific[:Fhist]) == length(LMTR_out.solver_specific[:NLSGradHist]) + @test LMTR_out.solver_specific[:NLSGradHist][end] == + bpdn_nls2.counters.neval_jprod_residual + bpdn_nls2.counters.neval_jtprod_residual - 1 + @test obj(bpdn_nls2, LMTR_out.solution) == LMTR_out.solver_specific[:Fhist][end] + @test h(LMTR_out.solution) == LMTR_out.solver_specific[:Hhist][end] + @test LMTR_out.status == :first_order + end +end + +# LM with R2DH as subsolver +for (h, h_name) ∈ ((NormL0(λ), "l0"), ) + @testset "bpdn-with-bounds-ls-LM-$(h_name)-R2DH" begin + x0 = zeros(bpdn_nls2.meta.nvar) + @test has_bounds(bpdn_nls2) + LM_out = LM( + bpdn_nls2, + h, + options, + x0 = x0, + subsolver = R2DH, + subsolver_options = subsolver_options, + ) + @test typeof(LM_out.solution) == typeof(bpdn_nls2.meta.x0) + @test length(LM_out.solution) == bpdn_nls2.meta.nvar + @test typeof(LM_out.solver_specific[:Fhist]) == typeof(LM_out.solution) + @test typeof(LM_out.solver_specific[:Hhist]) == typeof(LM_out.solution) + @test typeof(LM_out.solver_specific[:SubsolverCounter]) == Array{Int, 1} + @test typeof(LM_out.dual_feas) == eltype(LM_out.solution) + @test length(LM_out.solver_specific[:Fhist]) == length(LM_out.solver_specific[:Hhist]) + @test length(LM_out.solver_specific[:Fhist]) == + length(LM_out.solver_specific[:SubsolverCounter]) + @test length(LM_out.solver_specific[:Fhist]) == length(LM_out.solver_specific[:NLSGradHist]) + @test LM_out.solver_specific[:NLSGradHist][end] == + bpdn_nls2.counters.neval_jprod_residual + bpdn_nls2.counters.neval_jtprod_residual - 1 + @test obj(bpdn_nls2, LM_out.solution) == LM_out.solver_specific[:Fhist][end] + @test h(LM_out.solution) == LM_out.solver_specific[:Hhist][end] + @test LM_out.status == :first_order end end