Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Feb 26, 2025
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
2130fc9
update LM solver with possibility of using R2DH as subsolver and with…
MohamedLaghdafHABIBOULLAH Sep 27, 2024
825729d
Apply suggestions from code review
MohamedLaghdafHABIBOULLAH Feb 20, 2025
435a069
correct the definition of sqrt_ξ1_νInv
MohamedLaghdafHABIBOULLAH Feb 20, 2025
1787811
update νInv according to Maxence Proposal
MohamedLaghdafHABIBOULLAH Feb 20, 2025
b1c8037
remove theta update, as it will be rebased later
MohamedLaghdafHABIBOULLAH Feb 24, 2025
579cf06
declare sqrt_xi_1 as a local variable and add unit test for LM with R…
MohamedLaghdafHABIBOULLAH Feb 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
31 changes: 19 additions & 12 deletions src/LM_alg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.ν
Expand All @@ -85,7 +87,6 @@ function LM(
end

# initialize parameters
σk = max(1 / options.ν, σmin)
xk = copy(x0)
hk = h(xk[selected])
if hk == Inf
Expand All @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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])
Expand Down
88 changes: 63 additions & 25 deletions test/test_bounds.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading