LM: only recompute Jacobian norm if problem is nonlinear (#156)

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

Author: MohamedLaghdafHABIBOULLAH

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
@@ -101,6 +102,7 @@ function LM(
   xkn = similar(xk)
 
   local ξ1
+  local sqrt_ξ1_νInv
   k = 0
   Fobj_hist = zeros(maxIter)
   Hobj_hist = zeros(maxIter)
@@ -110,7 +112,7 @@ function LM(
 
   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 +179,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 = sqrt(ξ1 * νInv)
 
     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 +226,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 +248,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
@@ -262,9 +269,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 +288,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])

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