Skip to content

Add monotonicity parameter to LMSolver and TRSolver #226

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
dpo merged 3 commits into from
Sep 29, 2025
Merged

Add monotonicity parameter to LMSolver and TRSolver #226

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2ff471b
Add monotonicity parameter to LMSolver and TRSolver with history trac…
MohamedLaghdafHABIBOULLAH Sep 26, 2025
8b5abdc
Update src/LM_alg.jl
MohamedLaghdafHABIBOULLAH Sep 26, 2025
04beb5d
Update src/LM_alg.jl
MohamedLaghdafHABIBOULLAH Sep 27, 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
22 changes: 18 additions & 4 deletions src/LM_alg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,12 +24,13 @@ mutable struct LMSolver{
u_bound::V
l_bound_m_x::V
u_bound_m_x::V
m_fh_hist::V
subsolver::ST
subpb::PB
substats::GenericExecutionStats{T, V, V, T}
end

function LMSolver(reg_nls::AbstractRegularizedNLPModel{T, V}; subsolver = R2Solver) where {T, V}
function LMSolver(reg_nls::AbstractRegularizedNLPModel{T, V}; subsolver = R2Solver, m_monotone::Int = 1) where {T, V}
x0 = reg_nls.model.meta.x0
l_bound = reg_nls.model.meta.lvar
u_bound = reg_nls.model.meta.uvar
Expand All @@ -54,6 +55,8 @@ function LMSolver(reg_nls::AbstractRegularizedNLPModel{T, V}; subsolver = R2Solv
u_bound_m_x = similar(xk, 0)
end

m_fh_hist = fill(T(-Inf), m_monotone - 1)

ψ =
has_bnds ? shifted(reg_nls.h, xk, l_bound_m_x, u_bound_m_x, reg_nls.selected) :
shifted(reg_nls.h, xk)
Expand All @@ -80,6 +83,7 @@ function LMSolver(reg_nls::AbstractRegularizedNLPModel{T, V}; subsolver = R2Solv
u_bound,
l_bound_m_x,
u_bound_m_x,
m_fh_hist,
subsolver,
subpb,
substats,
Expand All @@ -105,7 +109,7 @@ where F(x) and J(x) are the residual and its Jacobian at x, respectively, ψ(s;

For advanced usage, first define a solver "LMSolver" to preallocate the memory used in the algorithm, and then call `solve!`:

solver = LMSolver(reg_nls; subsolver = R2Solver)
solver = LMSolver(reg_nls; subsolver = R2Solver, m_monotone = 1)
solve!(solver, reg_nls)

stats = RegularizedExecutionStats(reg_nls)
Expand All @@ -130,6 +134,7 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u
- `η2::T = T(0.9)`: very successful iteration threshold;
- `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful;
- `θ::T = 1/(1 + eps(T)^(1 / 5))`: is the model decrease fraction with respect to the decrease of the Cauchy model;
- `m_monotone::Int = 1`: monotonicity parameter. By default, LM is monotone but the non-monotone variant will be used if `m_monotone > 1`;
- `subsolver = R2Solver`: the solver used to solve the subproblems.

The algorithm stops either when `√(ξₖ/νₖ) < atol + rtol*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure.
Expand Down Expand Up @@ -166,7 +171,8 @@ end
function LM(reg_nls::AbstractRegularizedNLPModel; kwargs...)
kwargs_dict = Dict(kwargs...)
subsolver = pop!(kwargs_dict, :subsolver, R2Solver)
solver = LMSolver(reg_nls, subsolver = subsolver)
m_monotone = pop!(kwargs_dict, :m_monotone, 1)
solver = LMSolver(reg_nls, subsolver = subsolver, m_monotone = m_monotone)
stats = RegularizedExecutionStats(reg_nls)
solve!(solver, reg_nls, stats; kwargs_dict...)
return stats
Expand Down Expand Up @@ -214,7 +220,11 @@ function SolverCore.solve!(
ψ = solver.ψ
xkn = solver.xkn
s = solver.s
m_fh_hist = solver.m_fh_hist .= T(-Inf)
has_bnds = solver.has_bnds

m_monotone = length(m_fh_hist) + 1

if has_bnds
l_bound = solver.l_bound
u_bound = solver.u_bound
Expand Down Expand Up @@ -275,6 +285,7 @@ function SolverCore.solve!(
set_solver_specific!(stats, :smooth_obj, fk)
set_solver_specific!(stats, :nonsmooth_obj, hk)
set_solver_specific!(stats, :prox_evals, prox_evals + 1)
m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)

φ1 = let Fk = Fk, ∇fk = ∇fk
d -> dot(Fk, Fk) / 2 + dot(∇fk, d) # ∇fk = Jk^T Fk
Expand Down Expand Up @@ -345,7 +356,8 @@ function SolverCore.solve!(
error("LM: failed to compute a step: ξ = $ξ")
end

Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk
Δobj = fhmax - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
ρk = Δobj / ξ

verbose > 0 &&
Expand Down Expand Up @@ -400,6 +412,8 @@ function SolverCore.solve!(
σk = σk * γ
end

m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)

set_objective!(stats, fk + hk)
set_solver_specific!(stats, :smooth_obj, fk)
set_solver_specific!(stats, :nonsmooth_obj, hk)
Expand Down
20 changes: 17 additions & 3 deletions src/TR_alg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ mutable struct TRSolver{
u_bound::V
l_bound_m_x::V
u_bound_m_x::V
m_fh_hist::V
subsolver::ST
subpb::PB
substats::GenericExecutionStats{T, V, V, T}
Expand All @@ -32,6 +33,7 @@ function TRSolver(
reg_nlp::AbstractRegularizedNLPModel{T, V};
χ::X = NormLinf(one(T)),
subsolver = R2Solver,
m_monotone::Int = 1,
) where {T, V, X}
x0 = reg_nlp.model.meta.x0
l_bound = reg_nlp.model.meta.lvar
Expand All @@ -54,6 +56,8 @@ function TRSolver(
u_bound_m_x = similar(xk, 0)
end

m_fh_hist = fill(T(-Inf), m_monotone - 1)

ψ =
has_bnds || subsolver == TRDHSolver ?
shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) :
Expand Down Expand Up @@ -81,6 +85,7 @@ function TRSolver(
u_bound,
l_bound_m_x,
u_bound_m_x,
m_fh_hist,
subsolver,
subpb,
substats,
Expand All @@ -107,7 +112,7 @@ where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs + ½ sᵀ Bₖ s is a quadratic a

For advanced usage, first define a solver "TRSolver" to preallocate the memory used in the algorithm, and then call `solve!`:

solver = TR(reg_nlp; χ = NormLinf(1), subsolver = R2Solver)
solver = TRSolver(reg_nlp; χ = NormLinf(1), subsolver = R2Solver, m_monotone = 1)
solve!(solver, reg_nlp)

stats = RegularizedExecutionStats(reg_nlp)
Expand All @@ -129,6 +134,7 @@ For advanced usage, first define a solver "TRSolver" to preallocate the memory u
- `η1::T = √√eps(T)`: successful iteration threshold;
- `η2::T = T(0.9)`: very successful iteration threshold;
- `γ::T = T(3)`: trust-region radius parameter multiplier. Must satisfy `γ > 1`. The trust-region radius is updated as Δ := Δ*γ when the iteration is very successful and Δ := Δ/γ when the iteration is unsuccessful;
- `m_monotone::Int = 1`: monotonicity parameter. By default, TR is monotone but the non-monotone variant will be used if `m_monotone > 1`;
- `χ::F = NormLinf(1)`: norm used to define the trust-region;`
- `subsolver::S = R2Solver`: subsolver used to solve the subproblem that appears at each iteration.
- `sub_kwargs::NamedTuple = NamedTuple()`: a named tuple containing the keyword arguments to be sent to the subsolver. The solver will fail if invalid keyword arguments are provided to the subsolver. For example, if the subsolver is `R2Solver`, you can pass `sub_kwargs = (max_iter = 100, σmin = 1e-6,)`.
Expand Down Expand Up @@ -175,7 +181,8 @@ function TR(reg_nlp::AbstractRegularizedNLPModel{T, V}; kwargs...) where {T, V}
kwargs_dict = Dict(kwargs...)
subsolver = pop!(kwargs_dict, :subsolver, R2Solver)
χ = pop!(kwargs_dict, :χ, NormLinf(one(T)))
solver = TRSolver(reg_nlp, subsolver = subsolver, χ = χ)
m_monotone = pop!(kwargs_dict, :m_monotone, 1)
solver = TRSolver(reg_nlp, subsolver = subsolver, χ = χ, m_monotone = m_monotone)
stats = RegularizedExecutionStats(reg_nlp)
solve!(solver, reg_nlp, stats; kwargs_dict...)
return stats
Expand Down Expand Up @@ -221,8 +228,11 @@ function SolverCore.solve!(
xkn = solver.xkn
s = solver.s
χ = solver.χ
m_fh_hist = solver.m_fh_hist .= T(-Inf)
has_bnds = solver.has_bnds

m_monotone = length(m_fh_hist) + 1

if has_bnds || isa(solver.subsolver, TRDHSolver) #TODO elsewhere ?
l_bound_m_x, u_bound_m_x = solver.l_bound_m_x, solver.u_bound_m_x
l_bound, u_bound = solver.l_bound, solver.u_bound
Expand Down Expand Up @@ -293,6 +303,7 @@ function SolverCore.solve!(
set_solver_specific!(stats, :smooth_obj, fk)
set_solver_specific!(stats, :nonsmooth_obj, hk)
set_solver_specific!(stats, :prox_evals, prox_evals + 1)
m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)
Copy link

Copilot AI Sep 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[nitpick] The history update logic is duplicated in two places. Consider extracting this into a helper function or variable to reduce code duplication and improve maintainability.

Copilot uses AI. Check for mistakes.
Copy link
Collaborator Author

@MohamedLaghdafHABIBOULLAH MohamedLaghdafHABIBOULLAH Sep 26, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe in a future PR


# models
φ1 = let ∇fk = ∇fk
Expand Down Expand Up @@ -380,7 +391,8 @@ function SolverCore.solve!(
hkn = @views h(xkn[selected])
sNorm = χ(s)

Δobj = fk + hk - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk
Δobj = fhmax - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps()
ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps()

if (ξ ≤ 0 || isnan(ξ))
Expand Down Expand Up @@ -452,6 +464,8 @@ function SolverCore.solve!(
end
end

m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk)

set_objective!(stats, fk + hk)
set_solver_specific!(stats, :smooth_obj, fk)
set_solver_specific!(stats, :nonsmooth_obj, hk)
Expand Down
Loading