LMTR with bounds

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()

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,19 +84,21 @@ 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
   hk == -Inf && error("nonsmooth term is not proper")
 
   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

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