From 0c957d1d25b7761481833803dcd096712c86daff Mon Sep 17 00:00:00 2001 From: Geoffroy Leconte Date: Thu, 9 Feb 2023 16:20:19 -0500 Subject: [PATCH] LMTR with bounds --- examples/demo-bpdn-constr.jl | 12 +++++++++--- src/LMTR_alg.jl | 26 +++++++++++++++++++------- test/test_bounds.jl | 24 ++++++++++++++++++++++++ 3 files changed, 52 insertions(+), 10 deletions(-) diff --git a/examples/demo-bpdn-constr.jl b/examples/demo-bpdn-constr.jl index 2eb3a12c..4f72143e 100644 --- a/examples/demo-bpdn-constr.jl +++ b/examples/demo-bpdn-constr.jl @@ -8,7 +8,7 @@ include("plot-utils-bpdn.jl") Random.seed!(1234) -function demo_solver(f, sol, h, χ, suffix = "l0-linf") +function demo_solver(f, nls, sol, h, χ, suffix = "l0-linf") options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10) @info " using TR to solve with" h χ @@ -22,15 +22,21 @@ function demo_solver(f, sol, h, χ, suffix = "l0-linf") R2_out = R2(f, h, options, x0 = f.meta.x0) @info "R2 relative error" norm(R2_out.solution - sol) / norm(sol) plot_bpdn(R2_out, sol, "constr-r2-$(suffix)") + + @info " using LMTR to solve with" h χ + reset!(nls) + LMTR_out = LMTR(nls, h, χ, options, x0 = f.meta.x0) + @info "LMTR relative error" norm(LMTR_out.solution - sol) / norm(sol) + plot_bpdn(LMTR_out, sol, "constr-lmtr-r2-$(suffix)") end function demo_bpdn_constr(compound = 1) model, nls_model, sol = bpdn_model(compound, bounds = true) f = LSR1Model(model) λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10 - demo_solver(f, sol, NormL0(λ), NormLinf(1.0)) + demo_solver(f, nls_model, sol, NormL0(λ), NormLinf(1.0)) λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 3 - demo_solver(f, sol, NormL1(λ), NormLinf(1.0), "l1-linf") + demo_solver(f, nls_model, sol, NormL1(λ), NormLinf(1.0), "l1-linf") end demo_bpdn_constr() diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl index 1233ea30..b8425779 100644 --- a/src/LMTR_alg.jl +++ b/src/LMTR_alg.jl @@ -30,6 +30,7 @@ where F(x) and J(x) are the residual and its Jacobian at x, respectively, ψ(s; * `subsolver_logger::AbstractLogger`: a logger to pass to the subproblem solver * `subsolver`: the procedure used to compute a step (`PG` or `R2`) * `subsolver_options::ROSolverOptions`: default options to pass to the subsolver. +* `selected::AbstractVector{<:Integer}`: (default `1:f.meta.nvar`). ### Return values @@ -47,6 +48,7 @@ function LMTR( subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(), subsolver = R2, subsolver_options = ROSolverOptions(), + selected::AbstractVector{<:Integer} = 1:nls.meta.nvar, ) where {H, X} start_time = time() elapsed_time = 0.0 @@ -64,6 +66,12 @@ function LMTR( θ = options.θ β = options.β + local l_bound, u_bound + if has_bounds(nls) + l_bound = nls.meta.lvar + u_bound = nls.meta.uvar + end + if verbose == 0 ptf = Inf elseif verbose == 1 @@ -76,11 +84,11 @@ function LMTR( # initialize parameters xk = copy(x0) - hk = h(xk) + hk = h(xk[selected]) if hk == Inf verbose > 0 && @info "LMTR: finding initial guess where nonsmooth term is finite" prox!(xk, h, x0, one(eltype(x0))) - hk = h(xk) + hk = h(xk[selected]) hk < Inf || error("prox computation must be erroneous") verbose > 0 && @debug "LMTR: found point where h has value" hk end @@ -88,7 +96,9 @@ function LMTR( xkn = similar(xk) s = zero(xk) - ψ = shifted(h, xk, Δk, χ) + ψ = has_bounds(nls) ? shifted(h, xk, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk), selected) : + shifted(h, xk, Δk, χ) + Fobj_hist = zeros(maxIter) Hobj_hist = zeros(maxIter) @@ -168,7 +178,8 @@ function LMTR( subsolver_options.ϵa = k == 1 ? 1.0e-5 : max(ϵ, min(1.0e-1, ξ1 / 10)) ∆_effective = min(β * χ(s), Δk) - set_radius!(ψ, ∆_effective) + has_bounds(nls) ? set_bounds!(ψ, max.(-∆_effective, l_bound - xk), min.(∆_effective, u_bound - xk)) : + set_radius!(ψ, ∆_effective) s, iter, _ = with_logger(subsolver_logger) do subsolver(φ, ∇φ!, ψ, subsolver_options, s) end @@ -179,7 +190,7 @@ function LMTR( xkn .= xk .+ s residual!(nls, xkn, Fkn) fkn = dot(Fkn, Fkn) / 2 - hkn = h(xkn) + hkn = h(xkn[selected]) hkn == -Inf && error("nonsmooth term is not proper") Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps() @@ -203,11 +214,12 @@ function LMTR( if η2 ≤ ρk < Inf Δk = max(Δk, γ * sNorm) - set_radius!(ψ, Δk) + !has_bounds(nls) && set_radius!(ψ, Δk) end if η1 ≤ ρk < Inf xk .= xkn + has_bounds(nls) && set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk)) #update functions Fk .= Fkn @@ -224,7 +236,7 @@ function LMTR( if ρk < η1 || ρk == Inf Δk = Δk / 2 - set_radius!(ψ, Δk) + has_bounds(nls) ? set_bounds!(ψ, max.(-Δk, l_bound - xk), min.(Δk, u_bound - xk)) : set_radius!(ψ, Δk) end tired = k ≥ maxIter || elapsed_time > maxTime diff --git a/test/test_bounds.jl b/test/test_bounds.jl index 946f28c9..5cc3ca4e 100644 --- a/test/test_bounds.jl +++ b/test/test_bounds.jl @@ -25,3 +25,27 @@ for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "l end end end + +for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1")) + for solver_sym ∈ (:LMTR,) + solver_name = string(solver_sym) + solver = eval(solver_sym) + @testset "bpdn-with-bounds-ls-$(solver_name)-$(h_name)" 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) + @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