From 243e5279d4db4a099ab72cc29d4777bcfe921441 Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Thu, 30 Jul 2026 10:57:11 +0200 Subject: [PATCH 1/6] max-type nonmonotone R2 --- src/R2_alg.jl | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/src/R2_alg.jl b/src/R2_alg.jl index 96b99f34..8b39678a 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -18,6 +18,7 @@ mutable struct R2Solver{ u_bound::S l_bound_m_x::S u_bound_m_x::S + m_fh_hist::S Fobj_hist::Vector{R} Hobj_hist::Vector{R} Complex_hist::Vector{Int} @@ -29,6 +30,7 @@ function R2Solver( l_bound::S, u_bound::S; ψ = nothing, + m_monotone::Int = 6, ) where {R <: Real, S <: AbstractVector{R}} maxIter = options.maxIter xk = similar(x0) @@ -44,6 +46,7 @@ function R2Solver( l_bound_m_x = similar(xk, 0) u_bound_m_x = similar(xk, 0) end + m_fh_hist = fill(T(-Inf), m_monotone - 1) Fobj_hist = zeros(R, maxIter + 2) Hobj_hist = zeros(R, maxIter + 2) Complex_hist = zeros(Int, maxIter + 2) @@ -59,13 +62,18 @@ function R2Solver( u_bound, l_bound_m_x, u_bound_m_x, + m_fh_hist, Fobj_hist, Hobj_hist, Complex_hist, ) end -function R2Solver(reg_nlp::AbstractRegularizedNLPModel{T, V}; max_iter::Int = 10000) where {T, V} +function R2Solver( + reg_nlp::AbstractRegularizedNLPModel{T, V}; + max_iter::Int = 10000, + m_monotone::Int = 6, +) where {T, V} x0 = reg_nlp.model.meta.x0 l_bound = reg_nlp.model.meta.lvar u_bound = reg_nlp.model.meta.uvar @@ -85,6 +93,7 @@ function R2Solver(reg_nlp::AbstractRegularizedNLPModel{T, V}; max_iter::Int = 10 l_bound_m_x = similar(xk, 0) u_bound_m_x = similar(xk, 0) end + m_fh_hist = fill(T(-Inf), m_monotone - 1) Fobj_hist = zeros(T, max_iter + 2) Hobj_hist = zeros(T, max_iter + 2) Complex_hist = zeros(Int, max_iter + 2) @@ -104,6 +113,7 @@ function R2Solver(reg_nlp::AbstractRegularizedNLPModel{T, V}; max_iter::Int = 10 u_bound, l_bound_m_x, u_bound_m_x, + m_fh_hist, Fobj_hist, Hobj_hist, Complex_hist, @@ -129,7 +139,7 @@ where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs is the Taylor linear approximation For advanced usage, first define a solver "R2Solver" to preallocate the memory used in the algorithm, and then call `solve!`: - solver = R2Solver(reg_nlp) + solver = R2Solver(reg_nlp; m_monotone = 6) solve!(solver, reg_nlp) stats = RegularizedExecutionStats(reg_nlp) @@ -153,6 +163,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u - `η2::T = T(0.9)`: successful iteration threshold; - `ν::T = eps(T)^(1 / 5)`: multiplicative inverse of the regularization parameter: ν = 1/σ; - `γ::T = T(3)`: regularization parameter multiplier, σ := σ/γ when the iteration is very successful and σ := σγ when the iteration is unsuccessful. +- `m_monotone::Int = 6`: monotonicity parameter. By default, R2 is non-monotone but the monotone variant can be used with `m_monotone = 1`; - `compute_obj::Bool = true`: (advanced) whether `f(x₀)` should be computed or not. If set to false, then the value is retrieved from `stats.solver_specific[:smooth_obj]`; - `compute_grad::Bool = true`: (advanced) whether `∇f(x₀)` should be computed or not. If set to false, then the value is retrieved from `solver.∇fk`; @@ -288,8 +299,9 @@ end function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...) kwargs_dict = Dict(kwargs...) + m_monotone = pop!(kwargs_dict, :m_monotone, 6) max_iter = pop!(kwargs_dict, :max_iter, 10000) - solver = R2Solver(reg_nlp, max_iter = max_iter) + solver = R2Solver(reg_nlp, m_monotone = m_monotone, max_iter = max_iter) stats = GenericExecutionStats(reg_nlp.model) # TODO: change this to `stats = RegularizedExecutionStats(reg_nlp)` when FHist etc. is ruled out. cb = pop!( kwargs_dict, @@ -345,6 +357,7 @@ function SolverCore.solve!( ψ = solver.ψ xkn = solver.xkn s = solver.s + m_fh_hist = solver.m_fh_hist .= T(-Inf) has_bnds = solver.has_bnds if has_bnds l_bound, u_bound = solver.l_bound, solver.u_bound @@ -352,6 +365,7 @@ function SolverCore.solve!( update_bounds!(l_bound_m_x, u_bound_m_x, l_bound, u_bound, xk) set_bounds!(ψ, l_bound_m_x, u_bound_m_x) end + m_monotone = length(m_fh_hist) + 1 # initialize parameters improper = false @@ -401,6 +415,7 @@ function SolverCore.solve!( set_solver_specific!(stats, :smooth_obj, fk) set_solver_specific!(stats, :nonsmooth_obj, hk) set_solver_specific!(stats, :sigma, σk) + m_monotone > 1 && (m_fh_hist[stats.iter % (m_monotone - 1) + 1] = fk + hk) φk(d) = dot(∇fk, d) mk(d)::T = φk(d) + ψ(d)::T @@ -443,8 +458,10 @@ function SolverCore.solve!( hkn = @views h(xkn[selected]) improper = (hkn == -Inf) - Δobj = (fk + hk) - (fkn + hkn) + max(1, abs(fk + hk)) * 10 * eps() - ρk = Δobj / ξ + fhmax = m_monotone > 1 ? maximum(m_fh_hist) : fk + hk + Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps() + Δmod = fhmax - (fk + mks) + max(1, abs(hk)) * 10 * eps() + ρk = Δobj / Δmod verbose > 0 && stats.iter % verbose == 0 && @@ -486,6 +503,7 @@ function SolverCore.solve!( ν = 1 / σk @. mν∇fk = -ν * ∇fk + 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) From 3eef39342c5eff633af611314d24afb2aa54eec4 Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Thu, 30 Jul 2026 11:07:20 +0200 Subject: [PATCH 2/6] R2 without bookkeeping --- src/R2_alg.jl | 42 +++--------------------------------------- 1 file changed, 3 insertions(+), 39 deletions(-) diff --git a/src/R2_alg.jl b/src/R2_alg.jl index 8b39678a..b0961d49 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -19,9 +19,6 @@ mutable struct R2Solver{ l_bound_m_x::S u_bound_m_x::S m_fh_hist::S - Fobj_hist::Vector{R} - Hobj_hist::Vector{R} - Complex_hist::Vector{Int} end function R2Solver( @@ -32,7 +29,6 @@ function R2Solver( ψ = nothing, m_monotone::Int = 6, ) where {R <: Real, S <: AbstractVector{R}} - maxIter = options.maxIter xk = similar(x0) ∇fk = similar(x0) mν∇fk = similar(x0) @@ -47,9 +43,6 @@ function R2Solver( u_bound_m_x = similar(xk, 0) end m_fh_hist = fill(T(-Inf), m_monotone - 1) - Fobj_hist = zeros(R, maxIter + 2) - Hobj_hist = zeros(R, maxIter + 2) - Complex_hist = zeros(Int, maxIter + 2) return R2Solver( xk, ∇fk, @@ -63,15 +56,11 @@ function R2Solver( l_bound_m_x, u_bound_m_x, m_fh_hist, - Fobj_hist, - Hobj_hist, - Complex_hist, ) end function R2Solver( reg_nlp::AbstractRegularizedNLPModel{T, V}; - max_iter::Int = 10000, m_monotone::Int = 6, ) where {T, V} x0 = reg_nlp.model.meta.x0 @@ -94,9 +83,6 @@ function R2Solver( u_bound_m_x = similar(xk, 0) end m_fh_hist = fill(T(-Inf), m_monotone - 1) - Fobj_hist = zeros(T, max_iter + 2) - Hobj_hist = zeros(T, max_iter + 2) - Complex_hist = zeros(Int, max_iter + 2) ψ = has_bnds ? shifted(reg_nlp.h, xk, l_bound_m_x, u_bound_m_x, reg_nlp.selected) : @@ -114,9 +100,6 @@ function R2Solver( l_bound_m_x, u_bound_m_x, m_fh_hist, - Fobj_hist, - Hobj_hist, - Complex_hist, ) end @@ -241,9 +224,6 @@ function R2( kwargs..., ) outdict = Dict( - :Fhist => stats.solver_specific[:Fhist], - :Hhist => stats.solver_specific[:Hhist], - :Chist => stats.solver_specific[:SubsolverCounter], :NonSmooth => h, :status => stats.status, :fk => stats.solver_specific[:smooth_obj], @@ -284,9 +264,6 @@ function R2( kwargs..., ) outdict = Dict( - :Fhist => stats.solver_specific[:Fhist], - :Hhist => stats.solver_specific[:Hhist], - :Chist => stats.solver_specific[:SubsolverCounter], :NonSmooth => h, :status => stats.status, :fk => stats.solver_specific[:smooth_obj], @@ -300,22 +277,9 @@ end function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...) kwargs_dict = Dict(kwargs...) m_monotone = pop!(kwargs_dict, :m_monotone, 6) - max_iter = pop!(kwargs_dict, :max_iter, 10000) - solver = R2Solver(reg_nlp, m_monotone = m_monotone, max_iter = max_iter) - stats = GenericExecutionStats(reg_nlp.model) # TODO: change this to `stats = RegularizedExecutionStats(reg_nlp)` when FHist etc. is ruled out. - cb = pop!( - kwargs_dict, - :callback, - (nlp, solver, stats) -> begin - solver.Fobj_hist[stats.iter + 1] = stats.solver_specific[:smooth_obj] - solver.Hobj_hist[stats.iter + 1] = stats.solver_specific[:nonsmooth_obj] - solver.Complex_hist[stats.iter + 1] += 1 - end, - ) - solve!(solver, reg_nlp, stats; callback = cb, max_iter = max_iter, kwargs...) - set_solver_specific!(stats, :Fhist, solver.Fobj_hist[1:(stats.iter + 1)]) - set_solver_specific!(stats, :Hhist, solver.Hobj_hist[1:(stats.iter + 1)]) - set_solver_specific!(stats, :SubsolverCounter, solver.Complex_hist[1:(stats.iter + 1)]) + solver = R2Solver(reg_nlp, m_monotone = m_monotone) + stats = RegularizedExecutionStats(reg_nlp) + solve!(solver, reg_nlp, stats; kwargs...) return stats end From baa550efa4bc36a8b4069bae78485df6e5d34129 Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Thu, 30 Jul 2026 11:26:14 +0200 Subject: [PATCH 3/6] bugfix --- src/R2_alg.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/R2_alg.jl b/src/R2_alg.jl index b0961d49..00801f2c 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -279,7 +279,7 @@ function R2(reg_nlp::AbstractRegularizedNLPModel; kwargs...) m_monotone = pop!(kwargs_dict, :m_monotone, 6) solver = R2Solver(reg_nlp, m_monotone = m_monotone) stats = RegularizedExecutionStats(reg_nlp) - solve!(solver, reg_nlp, stats; kwargs...) + solve!(solver, reg_nlp, stats; kwargs_dict...) return stats end From 0d7d454227c77742b12f01bb7b68b2d69ad6012c Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Thu, 30 Jul 2026 11:27:32 +0200 Subject: [PATCH 4/6] test also monotone R2 --- test/runtests.jl | 6 ++++-- test/test_bounds.jl | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index caaf8e12..27f2cb8d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,11 +18,13 @@ const global bpdn, bpdn_nls, sol = bpdn_model(compound) const global bpdn2, bpdn_nls2, sol2 = bpdn_model(compound, bounds = true) const global λ = norm(grad(bpdn, zeros(bpdn.meta.nvar)), Inf) / 10 +R2mono(nlp, h, options; kwargs...) = R2(nlp, h, options; m_monotone=1, kwargs...) + include("test_AL.jl") for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs")) for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"), (IndBallL0(10 * compound), "B0")) - for solver_sym ∈ (:R2, :TR) + for solver_sym ∈ (:R2, :R2mono, :TR) solver_sym == :TR && mod_name == "exact" && continue solver_sym == :TR && h_name == "B0" && continue # FIXME solver_name = string(solver_sym) @@ -31,7 +33,7 @@ for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "l x0 = zeros(bpdn.meta.nvar) p = randperm(bpdn.meta.nvar)[1:nz] x0[p[1:nz]] = sign.(randn(nz)) # initial guess with nz nonzeros (necessary for h = B0) - args = solver_sym == :R2 ? () : (NormLinf(1.0),) + args = solver_sym ∈ (:R2, :R2mono) ? () : (NormLinf(1.0),) out = solver(mod(bpdn), h, args..., options, x0 = x0) @test typeof(out.solution) == typeof(bpdn.meta.x0) @test length(out.solution) == bpdn.meta.nvar diff --git a/test/test_bounds.jl b/test/test_bounds.jl index 0b173df8..dbe2ab8b 100644 --- a/test/test_bounds.jl +++ b/test/test_bounds.jl @@ -3,14 +3,14 @@ TR_TRDH(args...; kwargs...) = TR(args...; subsolver = TRDHSolver, kwargs...) for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs")) for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1")) - for solver_sym ∈ (:TR, :R2, :TR_TRDH) + for solver_sym ∈ (:TR, :R2, :R2mono, :TR_TRDH) solver_sym ∈ (:TR, :TR_TRDH) && mod_name == "exact" && continue solver_name = string(solver_sym) solver = eval(solver_sym) @testset "bpdn-with-bounds-$(mod_name)-$(solver_name)-$(h_name)" begin x0 = zeros(bpdn2.meta.nvar) p = randperm(bpdn2.meta.nvar)[1:nz] - args = solver_sym == :R2 ? () : (NormLinf(1.0),) + args = solver_sym ∈ (:R2, :R2mono) ? () : (NormLinf(1.0),) @test has_bounds(mod(bpdn2)) out = solver(mod(bpdn2), h, args..., options; x0 = x0) @test typeof(out.solution) == typeof(bpdn2.meta.x0) From a431e2b6a9744014f031e56d420708905e36ffd9 Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Fri, 31 Jul 2026 15:17:24 +0200 Subject: [PATCH 5/6] update bibtex --- docs/references.bib | 39 +++++++++++++++++++-------------------- 1 file changed, 19 insertions(+), 20 deletions(-) diff --git a/docs/references.bib b/docs/references.bib index 2baf1808..2204e0ba 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -31,26 +31,25 @@ @Article{leconte-orban-2025 doi = {10.1007/s10589-024-00604-5} } -@TechReport{diouane-gollier-orban-2024, - Author = {Diouane, Youssef and Gollier, Maxence and Orban, Dominique}, - Title = {A nonsmooth exact penalty method for equality-constrained optimization: complexity and implementation}, - Institution = {Groupe d’études et de recherche en analyse des décisions}, - Year = {2024}, - Type = {Les Cahiers du GERAD}, - Number = {G-2024-65}, - Address = {Montreal, Canada}, - doi = {10.48550/arxiv.2103.15993}, - url = {https://www.gerad.ca/fr/papers/G-2024-65}, +@Article{diouane-gollier-orban-2024, + author = {Diouane, Youssef and Gollier, Maxence and Orban, Dominique}, + title = {Nonsmooth exact penalty methods for equality-constrained optimization: complexity and implementation}, + journal = {SIAM Journal on Optimization}, + volume = {36}, + number = {2}, + pages = {626-650}, + year = {2026}, + doi = {10.1137/24M1705974} } -@TechReport{diouane-habiboullah-orban-2024, - Author = {Diouane, Youssef and Laghdaf Habiboullah, Mohamed and Orban, Dominique}, - Title = {A proximal modified quasi-Newton method for nonsmooth regularized optimization}, - Institution = {Groupe d’études et de recherche en analyse des décisions}, - Year = {2024}, - Type = {Les Cahiers du GERAD}, - Number = {G-2024-64}, - Address = {Montreal, Canada}, - doi = {10.48550/arxiv.2409.19428}, - url = {https://www.gerad.ca/fr/papers/G-2024-64}, +@Article{diouane-habiboullah-orban-2024, + author = {Diouane, Youssef and Laghdaf Habiboullah, Mohamed and Orban, Dominique}, + title = {A proximal modified quasi-Newton method for nonsmooth regularized optimization}, + journal = {SIAM Journal on Optimization}, + volume = {36}, + number = {2}, + pages = {534-563}, + year = {2026}, + doi = {10.1137/24M169761X} +} } \ No newline at end of file From 541be773740dcc016da9b3fdd7739d8cd018570f Mon Sep 17 00:00:00 2001 From: Alberto De Marchi Date: Fri, 31 Jul 2026 16:02:24 +0200 Subject: [PATCH 6/6] test R2 callback --- test/runtests.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 27f2cb8d..5bba396a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -142,5 +142,28 @@ for (mod, mod_name) ∈ ( end end +for (mod, mod_name) ∈ ((SpectralGradientModel, "spg"),) + for (h, h_name) ∈ ((NormL1(λ), "l1"),) + for solver_sym ∈ (:R2,:R2mono) + solver_name = string(solver_sym) + solver = eval(solver_sym) + @testset "bpdn-$(mod_name)-$(solver_name)-$(h_name)-callback" begin + x0 = zeros(bpdn.meta.nvar) + cb = (regnlp, solver, stats) -> begin + if stats.iter == 2 + stats.status = :user + end + end + out = solver(mod(bpdn), h, options, x0 = x0, callback=cb) + @test typeof(out.solution) == typeof(bpdn.meta.x0) + @test length(out.solution) == bpdn.meta.nvar + @test typeof(out.dual_feas) == eltype(out.solution) + @test out.status == :user + @test out.iter == 2 + end + end + end +end + include("test_bounds.jl") include("test_allocs.jl")