From ea6fa12a90d4ee74c7fee6db2f75510c40eb392b Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 11:08:33 -0500 Subject: [PATCH 1/9] R2N & R2DH: add cauchy step stopping criterion --- src/R2DH.jl | 41 +++++++++++++++++++++++++++++++---------- src/R2N.jl | 40 ++++++++++++++++++++++++++++++---------- test/runtests.jl | 14 ++++++++++++++ 3 files changed, 75 insertions(+), 20 deletions(-) diff --git a/src/R2DH.jl b/src/R2DH.jl index de93ffcb..507761c1 100644 --- a/src/R2DH.jl +++ b/src/R2DH.jl @@ -121,6 +121,10 @@ or - `x::V = nlp.meta.x0`: the initial guess; - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -136,7 +140,7 @@ or - `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`; -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. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ` is the current step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -168,7 +172,7 @@ function R2DH( η1 = options.η1, η2 = options.η2, γ = options.γ, - θ = options.θ, + θ = options.θ; kwargs_dict..., ) end @@ -225,6 +229,10 @@ function SolverCore.solve!( x::V = reg_nlp.model.meta.x0, atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = eps(T)^(1 / 4), verbose::Int = 0, max_iter::Int = 10000, @@ -286,12 +294,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:iter, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :arrow], - [Int, T, T, T, T, T, T, T, Char], + [:iter, :fx, :hx, :xi, :normsdnu, :ρ, :σ, :normx, :norms, :arrow], + [Int, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( # TODO: Add this as constant dict elsewhere :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ/ν)", + :normsdnu => "‖sₖ‖/ν", :normx => "‖x‖", :norms => "‖s‖", :arrow => "R2DH", @@ -341,12 +350,18 @@ function SolverCore.solve!( mks = mk(s) + # Estimate optimality and check stopping criteria ξ = hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν₁) : sqrt(-ξ / ν₁) - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2DH: prox-gradient step should produce a decrease but ξ = $(ξ)") - atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ_νInv # make stopping test absolute and relative + + norm_s = norm(s) + norm_sdν = norm_s / ν₁ + atol_step += rtol_step * norm_sdν # make stopping test absolute and relative + + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) set_status!( stats, @@ -386,10 +401,11 @@ function SolverCore.solve!( fk, hk, sqrt_ξ_νInv, + norm_sdν, ρk, σk, norm(xk), - norm(s), + norm_s, (η2 ≤ ρk < Inf) ? '↘' : (ρk < η1 ? '↗' : '='), ], colsep = 1, @@ -439,11 +455,16 @@ function SolverCore.solve!( spectral_test ? prox!(s, ψ, mν∇fk, ν₁) : iprox!(s, ψ, ∇fk, dkσk) mks = mk(s) + # Estimate optimality and check stopping criteria ξ = hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν₁) : sqrt(-ξ / ν₁) - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2DH: prox-gradient step should produce a decrease but ξ = $(ξ)") + + norm_s = norm(s) + norm_sdν = norm_s / ν₁ + + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) set_status!( stats, @@ -465,8 +486,8 @@ function SolverCore.solve!( end if verbose > 0 && stats.status == :first_order - @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), ""], colsep = 1) - @info "R2DH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)" + @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, norm_sdν, ρk, σk, norm(xk), norm_s, ""], colsep = 1) + @info "R2DH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖sₖ‖/ν = $(norm_sdν)" end set_solution!(stats, xk) diff --git a/src/R2N.jl b/src/R2N.jl index 51ad3c56..e04a521b 100644 --- a/src/R2N.jl +++ b/src/R2N.jl @@ -136,6 +136,10 @@ For advanced usage, first define a solver "R2NSolver" to preallocate the memory - `x::V = nlp.meta.x0`: the initial guess; - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance; - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -153,7 +157,7 @@ For advanced usage, first define a solver "R2NSolver" to preallocate the memory - `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`; - `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,)`. -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. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -215,6 +219,10 @@ function SolverCore.solve!( x::V = reg_nlp.model.meta.x0, atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = eps(T)^(1 / 4), verbose::Int = 0, max_iter::Int = 10000, @@ -278,12 +286,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:outer, :inner, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :normB, :arrow], - [Int, Int, T, T, T, T, T, T, T, T, Char], + [:outer, :inner, :fx, :hx, :xi, :norms1dnu, :ρ, :σ, :normx, :norms, :normB, :arrow], + [Int, Int, T, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ1/ν)", + :norms1dnu => "‖sₖ₁‖/ν", :normx => "‖x‖", :norms => "‖s‖", :normB => "‖B‖", @@ -344,12 +353,18 @@ function SolverCore.solve!( prox!(s1, ψ, mν∇fk, ν₁) mks = mk1(s1) + # Estimate optimality and check stopping criteria ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("R2N: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - atol += rtol * sqrt_ξ1_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ1_νInv # make stopping test absolute and relative + + norm_s_cauchy = norm(s1) + norm_s_cauchydν = norm_s_cauchy / ν₁ + atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) set_status!( stats, @@ -428,6 +443,7 @@ function SolverCore.solve!( fk, hk, sqrt_ξ1_νInv, + norm_s_cauchydν, ρk, σk, norm(xk), @@ -491,13 +507,17 @@ function SolverCore.solve!( prox!(s1, ψ, mν∇fk, ν₁) mks = mk1(s1) + # Estimate optimality and check stopping criteria ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps() - sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) - (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("R2N: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") + + norm_s_cauchy = norm(s1) + norm_s_cauchydν = norm_s_cauchy / ν₁ + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) + set_status!( stats, get_status( @@ -519,10 +539,10 @@ function SolverCore.solve!( if verbose > 0 && stats.status == :first_order @info log_row( - Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, ρk, σk, norm(xk), norm(s), λmax, ""], + Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, σk, norm(xk), norm(s), λmax, ""], colsep = 1, ) - @info "R2N: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" + @info "R2N: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖sₖ₁‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) diff --git a/test/runtests.jl b/test/runtests.jl index caaf8e12..3d0a6ee2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -135,6 +135,20 @@ for (mod, mod_name) ∈ ( @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order @test out.step_status == (out.iter > 0 ? :accepted : :unknown) + + # Test with the different stopping criteria + out = solver(mod(bpdn), h, options; x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + @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 == :first_order + + out = solver(mod(bpdn), h, options; x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + @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 == :first_order + end end end From 2f7381672973772c0a77f3565d8db9b1dc1b6d46 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 12:12:12 -0500 Subject: [PATCH 2/9] TR & TRDH: add cauchy step stopping criterion --- src/TRDH_alg.jl | 33 ++++++++++++++++++++++++--------- src/TR_alg.jl | 46 +++++++++++++++++++++++++++++++++------------- test/runtests.jl | 30 ++++++++++++++++++++++++++++-- 3 files changed, 85 insertions(+), 24 deletions(-) diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl index 6a76ce0c..ef50c6a9 100644 --- a/src/TRDH_alg.jl +++ b/src/TRDH_alg.jl @@ -131,6 +131,10 @@ For advanced usage, first define a solver "TRDHSolver" to preallocate the memory - `x::V = nlp.meta.x0`: the initial guess; - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance; - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -146,7 +150,7 @@ For advanced usage, first define a solver "TRDHSolver" to preallocate the memory - `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`; -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. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ` is the current step. Alternatively, if `reduce_TR = true`, then ξₖ₁ := f(xₖ) + h(xₖ) - φ(sₖ₁; xₖ) - ψ(sₖ₁; xₖ) is used instead of ξₖ, where sₖ₁ is the Cauchy point. # Output @@ -239,6 +243,10 @@ function SolverCore.solve!( x::V = reg_nlp.model.meta.x0, atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = eps(T)^(1 / 4), verbose::Int = 0, max_iter::Int = 10000, @@ -372,10 +380,9 @@ function SolverCore.solve!( ξ1 = hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol) (ξ1 < 0 && sqrt_ξ_νInv > neg_tol) && error("TR: prox-gradient step should produce a decrease but ξ = $(ξ)") - atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ_νInv # make stopping test absolute and relative end Δ_effective = reduce_TR ? min(β * χ(s), Δk) : Δk @@ -394,12 +401,16 @@ function SolverCore.solve!( if !reduce_TR sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("TRDH: prox-gradient step should produce a decrease but ξ = $(ξ)") - atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative #TODO : this is redundant code with the other case of the test. + atol_decr += rtol_decr * sqrt_ξ_νInv # make stopping test absolute and relative #TODO : this is redundant code with the other case of the test. end + norm_s = norm(s) + norm_sdν = norm_s / ν + atol_step += rtol_step * norm_sdν # make stopping test absolute and relative + + solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) | (norm_sdν ≤ atol_step) set_status!( stats, get_status( @@ -434,6 +445,7 @@ function SolverCore.solve!( fk, hk, sqrt_ξ_νInv, + norm_sdν, ρk, Δk, χ(xk), @@ -491,7 +503,6 @@ function SolverCore.solve!( prox!(s, ψ, mν∇fk, ν) ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv < atol) (ξ1 < 0 && sqrt_ξ_νInv > neg_tol) && error("TRDH: prox-gradient step should produce a decrease but ξ = $(ξ)") end @@ -503,11 +514,15 @@ function SolverCore.solve!( if !reduce_TR ξ = hk - mk(s) + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("TRDH: prox-gradient step should produce a decrease but ξ = $(ξ)") end + norm_s = norm(s) + norm_sdν = norm_s / ν + + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) + set_status!( stats, get_status( @@ -529,10 +544,10 @@ function SolverCore.solve!( if verbose > 0 && stats.status == :first_order @info log_row( - Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, Δk, χ(xk), sNorm, norm(D.d), ""], + Any[stats.iter, fk, hk, sqrt_ξ_νInv, norm_sdν, ρk, Δk, χ(xk), sNorm, norm(D.d), ""], colsep = 1, ) - @info "TRDH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)" + @info "TRDH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖sₖ‖/ν = $(norm_sdν)" end set_solution!(stats, xk) diff --git a/src/TR_alg.jl b/src/TR_alg.jl index 7d3d9bbe..8c1e6dbc 100644 --- a/src/TR_alg.jl +++ b/src/TR_alg.jl @@ -136,6 +136,10 @@ For advanced usage, first define a solver "TRSolver" to preallocate the memory u - `x::V = nlp.meta.x0`: the initial guess; - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance (see stopping conditions below); - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -153,7 +157,7 @@ For advanced usage, first define a solver "TRSolver" to preallocate the memory u - `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`; - `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,)`. -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. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -211,6 +215,10 @@ function SolverCore.solve!( atol::T = √eps(T), sub_atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = eps(T)^(1 / 4), verbose::Int = 0, subsolver_logger::Logging.AbstractLogger = Logging.SimpleLogger(), @@ -276,12 +284,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:outer, :inner, :fx, :hx, :xi, :ρ, :Δ, :normx, :norms, :normB, :arrow], - [Int, Int, T, T, T, T, T, T, T, T, Char], + [:outer, :inner, :fx, :hx, :xi, :norms1dnu, :ρ, :Δ, :normx, :norms, :normB, :arrow], + [Int, Int, T, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( # TODO: Add this as constant dict elsewhere :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ1/ν)", + :norms1dnu => "‖sₖ₁‖/ν", :normx => "‖x‖", :norms => "‖s‖", :normB => "‖B‖", @@ -340,16 +349,21 @@ function SolverCore.solve!( end prox!(s, ψ, mν∇fk, ν₁) - ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() - ξ1 > 0 || error("TR: first prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - sqrt_ξ1_νInv = sqrt(ξ1 / ν₁) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) + # Estimate optimality and check stopping criteria + ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() + sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("TR: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - atol += rtol * sqrt_ξ1_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ1_νInv # make stopping test absolute and relative sub_atol += rtol * sqrt_ξ1_νInv + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν₁ + atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) + set_status!( stats, get_status( @@ -431,6 +445,7 @@ function SolverCore.solve!( fk, hk, sqrt_ξ1_νInv, + norm_s_cauchydν, ρk, Δk, χ(xk), @@ -504,12 +519,17 @@ function SolverCore.solve!( @. mν∇fk = -ν₁ * ∇fk prox!(s, ψ, mν∇fk, ν₁) - ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() - sqrt_ξ1_νInv = sqrt(ξ1 / ν₁) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) + # Estimate optimality and check stopping criteria + ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() + sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν₁) : sqrt(-ξ1 / ν₁) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("TR: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") + + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν₁ + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) set_status!( stats, @@ -531,10 +551,10 @@ function SolverCore.solve!( end if verbose > 0 && stats.status == :first_order @info log_row( - Any[stats.iter, solver.substats.iter, fk, hk, sqrt_ξ1_νInv, ρk, Δk, χ(xk), χ(s), λmax, ""], + Any[stats.iter, solver.substats.iter, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, Δk, χ(xk), χ(s), λmax, ""], colsep = 1, ) - @info "TR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" + @info "TR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖sₖ₁‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) diff --git a/test/runtests.jl b/test/runtests.jl index 3d0a6ee2..4f34c45c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -57,6 +57,19 @@ for (mod, mod_name) ∈ ((SpectralGradientModel, "spg"),) @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order @test out.step_status == (out.iter > 0 ? :accepted : :unknown) + + # Test with the different stopping criteria + out = TRDH(mod(bpdn), h, χ, options, x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + @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 == :first_order + + out = TRDH(mod(bpdn), h, χ, options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + @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 == :first_order end end end @@ -74,6 +87,19 @@ for (mod, mod_name) ∈ ((LSR1Model, "lsr1"), (LBFGSModel, "lbfgs")) @test typeof(TR_out.dual_feas) == eltype(TR_out.solution) @test TR_out.status == :first_order @test TR_out.step_status == (TR_out.iter > 0 ? :accepted : :unknown) + + # Test with the different stopping criteria + TR_out = TR(mod(bpdn), h, NormL2(1.0), options, x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + @test typeof(TR_out.solution) == typeof(bpdn.meta.x0) + @test length(TR_out.solution) == bpdn.meta.nvar + @test typeof(TR_out.dual_feas) == eltype(TR_out.solution) + @test TR_out.status == :first_order + + TR_out = TR(mod(bpdn), h, NormL2(1.0), options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + @test typeof(TR_out.solution) == typeof(bpdn.meta.x0) + @test length(TR_out.solution) == bpdn.meta.nvar + @test typeof(TR_out.dual_feas) == eltype(TR_out.solution) + @test TR_out.status == :first_order end end end @@ -137,13 +163,13 @@ for (mod, mod_name) ∈ ( @test out.step_status == (out.iter > 0 ? :accepted : :unknown) # Test with the different stopping criteria - out = solver(mod(bpdn), h, options; x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + out = solver(mod(bpdn), h, options, x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) @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 == :first_order - out = solver(mod(bpdn), h, options; x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + out = solver(mod(bpdn), h, options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @test typeof(out.solution) == typeof(bpdn.meta.x0) @test length(out.solution) == bpdn.meta.nvar @test typeof(out.dual_feas) == eltype(out.solution) From 67ede8a99ce8a9c8e5e64e2b7561d49a1789c182 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 12:25:27 -0500 Subject: [PATCH 3/9] R2: add cauchy step stopping criterion --- src/R2_alg.jl | 40 +++++++++++++++++++++++++++++----------- test/runtests.jl | 13 +++++++++++++ 2 files changed, 42 insertions(+), 11 deletions(-) diff --git a/src/R2_alg.jl b/src/R2_alg.jl index 48bfc54e..fce199a7 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -143,6 +143,10 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u - `x::V = nlp.meta.x0`: the initial guess; - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ‖/ν₁` (see below); - `neg_tol::T = eps(T)^(1 / 4)`: negative tolerance - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -156,7 +160,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u - `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`; -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. +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 or when `‖sₖ‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ` is the current step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -315,6 +319,10 @@ function SolverCore.solve!( x::V = reg_nlp.model.meta.x0, atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = eps(T)^(1 / 4), verbose::Int = 0, max_iter::Int = 10000, @@ -367,13 +375,14 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:iter, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :arrow], - [Int, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Char], + [:iter, :fx, :hx, :xi, :normsdnu, :ρ, :σ, :normx, :norms, :arrow], + [Int, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( # TODO: Add this as constant dict elsewhere :iter => "iter", :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ/ν)", + :normsdnu => "‖s‖/ν", :ρ => "ρ", :σ => "σ", :normx => "‖x‖", @@ -408,15 +417,19 @@ function SolverCore.solve!( prox!(s, ψ, mν∇fk, ν) mks = mk(s) + # Estimate optimality and check stopping criteria ξ = hk - mks + max(1, abs(hk)) * 10 * eps() - sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) - atol += rtol * sqrt_ξ_νInv # make stopping test absolute and relative - - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)") + atol_decr += rtol_decr * sqrt_ξ_νInv # make stopping test absolute and relative + norm_s = norm(s) + norm_sdν = norm_s / ν + atol_step += rtol_step * norm_sdν # make stopping test absolute and relative + + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) + set_status!( stats, get_status( @@ -454,10 +467,11 @@ function SolverCore.solve!( fk, hk, sqrt_ξ_νInv, + norm_sdν, ρk, σk, norm(xk), - norm(s), + norm_s, (η2 ≤ ρk < Inf) ? '↘' : (ρk < η1 ? '↗' : '='), ], colsep = 1, @@ -499,10 +513,14 @@ function SolverCore.solve!( ξ = hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) - solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2: prox-gradient step should produce a decrease but ξ = $(ξ)") + norm_s = norm(s) + norm_sdν = norm_s / ν + + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) + set_status!( stats, get_status( @@ -523,8 +541,8 @@ function SolverCore.solve!( end if verbose > 0 && stats.status == :first_order - @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, ρk, σk, norm(xk), norm(s), ""], colsep = 1) - @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv)" + @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, norm_sdν, ρk, σk, norm(xk), norm_s, ""], colsep = 1) + @info "R2: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖s‖/ν = $(norm_sdν)" end set_solution!(stats, xk) diff --git a/test/runtests.jl b/test/runtests.jl index 4f34c45c..db3248c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -38,6 +38,19 @@ for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "l @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order @test out.step_status == (out.iter > 0 ? :accepted : :unknown) + + # Test with the different stopping criteria + out = solver(mod(bpdn), h, args..., options, x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + @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 == :first_order + + out = solver(mod(bpdn), h, args..., options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + @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 == :first_order end end end From 6d12c28794a511d9379320133dc3ff0fdf6fca43 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 12:39:29 -0500 Subject: [PATCH 4/9] fix logging --- src/R2DH.jl | 4 ++-- src/R2N.jl | 4 ++-- src/TRDH_alg.jl | 7 ++++--- src/TR_alg.jl | 2 +- 4 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/R2DH.jl b/src/R2DH.jl index 507761c1..7310832c 100644 --- a/src/R2DH.jl +++ b/src/R2DH.jl @@ -300,7 +300,7 @@ function SolverCore.solve!( :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ/ν)", - :normsdnu => "‖sₖ‖/ν", + :normsdnu => "‖s‖/ν", :normx => "‖x‖", :norms => "‖s‖", :arrow => "R2DH", @@ -487,7 +487,7 @@ function SolverCore.solve!( if verbose > 0 && stats.status == :first_order @info log_row(Any[stats.iter, fk, hk, sqrt_ξ_νInv, norm_sdν, ρk, σk, norm(xk), norm_s, ""], colsep = 1) - @info "R2DH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖sₖ‖/ν = $(norm_sdν)" + @info "R2DH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖s‖/ν = $(norm_sdν)" end set_solution!(stats, xk) diff --git a/src/R2N.jl b/src/R2N.jl index e04a521b..96b607f2 100644 --- a/src/R2N.jl +++ b/src/R2N.jl @@ -292,7 +292,7 @@ function SolverCore.solve!( :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ1/ν)", - :norms1dnu => "‖sₖ₁‖/ν", + :norms1dnu => "‖s₁‖/ν", :normx => "‖x‖", :norms => "‖s‖", :normB => "‖B‖", @@ -542,7 +542,7 @@ function SolverCore.solve!( Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, σk, norm(xk), norm(s), λmax, ""], colsep = 1, ) - @info "R2N: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖sₖ₁‖/ν = $(norm_s_cauchydν)" + @info "R2N: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖s₁‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl index ef50c6a9..2bb0362c 100644 --- a/src/TRDH_alg.jl +++ b/src/TRDH_alg.jl @@ -312,12 +312,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:iter, :fx, :hx, :xi, :ρ, :Δ, :normx, :norms, :normD, :arrow], - [Int, T, T, T, T, T, T, T, T, Char], + [:iter, :fx, :hx, :xi, :normsdnu, :ρ, :Δ, :normx, :norms, :normD, :arrow], + [Int, T, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( # TODO: Add this as constant dict elsewhere :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ/ν)", + :normsdnu => "‖s‖/ν", :normx => "‖x‖", :norms => "‖s‖", :normD => "‖D‖", @@ -547,7 +548,7 @@ function SolverCore.solve!( Any[stats.iter, fk, hk, sqrt_ξ_νInv, norm_sdν, ρk, Δk, χ(xk), sNorm, norm(D.d), ""], colsep = 1, ) - @info "TRDH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖sₖ‖/ν = $(norm_sdν)" + @info "TRDH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv) and ‖s‖/ν = $(norm_sdν)" end set_solution!(stats, xk) diff --git a/src/TR_alg.jl b/src/TR_alg.jl index 8c1e6dbc..8d7a3172 100644 --- a/src/TR_alg.jl +++ b/src/TR_alg.jl @@ -554,7 +554,7 @@ function SolverCore.solve!( Any[stats.iter, solver.substats.iter, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, Δk, χ(xk), χ(s), λmax, ""], colsep = 1, ) - @info "TR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖sₖ₁‖/ν = $(norm_s_cauchydν)" + @info "TR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖s₁‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) From f90adbc679c8d9c223938f39f0e6c7aa0e0574e9 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 12:52:29 -0500 Subject: [PATCH 5/9] LM & LMTR: add cauchy step stopping criterion --- src/LMTR_alg.jl | 38 ++++++++++++++++++++++++++++++-------- src/LM_alg.jl | 40 ++++++++++++++++++++++++++++++---------- test/runtests.jl | 13 +++++++++++++ 3 files changed, 73 insertions(+), 18 deletions(-) diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl index 26187432..afe70181 100644 --- a/src/LMTR_alg.jl +++ b/src/LMTR_alg.jl @@ -127,6 +127,10 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `atol::T = √eps(T)`: absolute tolerance; - `sub_atol::T = atol`: subsolver absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); - `neg_tol::T = zero(T): negative tolerance; - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -142,7 +146,7 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `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,)`. -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. +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 or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -193,6 +197,10 @@ function SolverCore.solve!( atol::T = √eps(T), sub_atol::T = atol, rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = zero(T), verbose::Int = 0, max_iter::Int = 10000, @@ -251,12 +259,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:outer, :inner, :fx, :hx, :xi, :ρ, :Δ, :normx, :norms, :ν, :arrow], - [Int, Int, T, T, T, T, T, T, T, T, Char], + [:outer, :inner, :fx, :hx, :xi, :norms1dnu, :ρ, :Δ, :normx, :norms, :ν, :arrow], + [Int, Int, T, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ1/ν)", + :norms1dnu => "‖s₁‖/ν", :normx => "‖x‖", :norms => "‖s‖", :arrow => "LMTR", @@ -302,14 +311,21 @@ function SolverCore.solve!( # Take first proximal gradient step s1 and see if current xk is nearly stationary. # s1 minimizes φ1(d) + ‖d‖² / 2 / ν + ψ(d) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0)) prox!(s, ψ, mν∇fk, ν) + + # Estimate optimality and check stopping criteria ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps() sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("LMTR: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - atol += rtol * sqrt_ξ1_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ1_νInv # make stopping test absolute and relative sub_atol += rtol * sqrt_ξ1_νInv + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν + atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) + set_status!( stats, get_status( @@ -390,6 +406,7 @@ function SolverCore.solve!( fk, hk, sqrt_ξ1_νInv, + norm_s_cauchydν, ρk, ∆_effective, χ(xk), @@ -455,12 +472,17 @@ function SolverCore.solve!( prox!(s, ψ, mν∇fk, ν) mks = mk1(s) + # Estimate optimality and check stopping criteria ξ1 = fk + hk - mks + max(1, abs(hk)) * 10 * eps() sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("LM: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) + set_status!( stats, get_status( @@ -481,8 +503,8 @@ function SolverCore.solve!( end if verbose > 0 && stats.status == :first_order - @info log_row(Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, ρk, Δk, χ(xk), χ(s), ν, ""], colsep = 1) - @info "LMTR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" + @info log_row(Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, Δk, χ(xk), χ(s), ν, ""], colsep = 1) + @info "LMTR: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖s₁‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) diff --git a/src/LM_alg.jl b/src/LM_alg.jl index ef514b44..93e0f6cd 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -127,6 +127,10 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `nonlinear::Bool = true`: whether the function `F` is nonlinear or not. - `atol::T = √eps(T)`: absolute tolerance; - `rtol::T = √eps(T)`: relative tolerance; +- `atol_decr::T = atol`: (advanced) absolute tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `rtol_decr::T = rtol`: (advanced) relative tolerance for the optimality measure `√(ξₖ/νₖ)` (see below); +- `atol_step::T = atol`: (advanced) absolute tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); +- `rtol_step::T = rtol`: (advanced) relative tolerance for the optimality measure `‖sₖ₁‖/ν₁` (see below); - `neg_tol::T = zero(T): negative tolerance; - `max_eval::Int = -1`: maximum number of evaluation of the objective function (negative number means unlimited); - `max_time::Float64 = 30.0`: maximum time limit in seconds; @@ -142,7 +146,7 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `subsolver = R2Solver`: the solver used to solve the subproblems. - `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,)`. -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. +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 or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀₁‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -192,6 +196,10 @@ function SolverCore.solve!( nonlinear::Bool = true, atol::T = √eps(T), rtol::T = √eps(T), + atol_decr::T = atol, + rtol_decr::T = rtol, + atol_step::T = atol, + rtol_step::T = rtol, neg_tol::T = zero(T), verbose::Int = 0, max_iter::Int = 10000, @@ -254,12 +262,13 @@ function SolverCore.solve!( if verbose > 0 @info log_header( - [:outer, :inner, :fx, :hx, :xi, :ρ, :σ, :normx, :norms, :normJ, :arrow], - [Int, Int, T, T, T, T, T, T, T, T, Char], + [:outer, :inner, :fx, :hx, :xi, :normsdnu, :ρ, :σ, :normx, :norms, :normJ, :arrow], + [Int, Int, T, T, T, T, T, T, T, T, T, Char], hdr_override = Dict{Symbol, String}( :fx => "f(x)", :hx => "h(x)", :xi => "√(ξ1/ν)", + :normsdnu => "‖s₁‖/ν", :normx => "‖x‖", :norms => "‖s‖", :normJ => "‖J‖²", @@ -306,12 +315,19 @@ function SolverCore.solve!( end prox!(s, ψ, mν∇fk, ν) + + # Estimate optimality and check stopping criteria ξ1 = fk + hk - mk1(s) + max(1, abs(fk + hk)) * 10 * eps() sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("LM: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") - atol += rtol * sqrt_ξ1_νInv # make stopping test absolute and relative + atol_decr += rtol_decr * sqrt_ξ1_νInv # make stopping test absolute and relative + + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν + atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) set_status!( stats, @@ -384,6 +400,7 @@ function SolverCore.solve!( fk, hk, sqrt_ξ1_νInv, + norm_s_cauchydν, ρk, σk, norm(xk), @@ -443,14 +460,17 @@ function SolverCore.solve!( prox!(s, ψ, mν∇fk, ν) mks = mk1(s) + # Estimate optimality and check stopping criteria ξ1 = fk + hk - mks + max(1, abs(hk)) * 10 * eps() - sqrt_ξ1_νInv = ξ1 ≥ 0 ? sqrt(ξ1 / ν) : sqrt(-ξ1 / ν) - solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol) - (ξ1 < 0 && sqrt_ξ1_νInv > neg_tol) && error("LM: prox-gradient step should produce a decrease but ξ1 = $(ξ1)") + norm_s_cauchy = norm(s) + norm_s_cauchydν = norm_s_cauchy / ν + + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) + set_status!( stats, get_status( @@ -472,10 +492,10 @@ function SolverCore.solve!( if verbose > 0 && stats.status == :first_order @info log_row( - Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, ρk, σk, norm(xk), norm(s), 1 / ν, ""], + Any[stats.iter, 0, fk, hk, sqrt_ξ1_νInv, norm_s_cauchydν, ρk, σk, norm(xk), norm(s), 1 / ν, ""], colsep = 1, ) - @info "LM: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" + @info "LM: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv) and ‖s‖/ν = $(norm_s_cauchydν)" end set_solution!(stats, xk) diff --git a/test/runtests.jl b/test/runtests.jl index db3248c8..2ceecf66 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -133,6 +133,19 @@ for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"), (IndBallL0(10 * com @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order @test out.step_status == (out.iter > 0 ? :accepted : :unknown) + + # Test with the different stopping criteria + out = solver(bpdn_nls, h, args..., options, x0 = x0, atol_decr = 1e-6, rtol_decr = 1e-6, atol_step = 0.0, rtol_step = 0.0) + @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 == :first_order + + out = solver(bpdn_nls, h, args..., options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) + @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 == :first_order end end end From ce0b8137b48b44349086e2f3c3bb59e607ed424c Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Fri, 20 Feb 2026 13:04:44 -0500 Subject: [PATCH 6/9] update set_residuals --- src/LMTR_alg.jl | 2 +- src/LM_alg.jl | 2 +- src/R2DH.jl | 2 +- src/R2N.jl | 2 +- src/R2_alg.jl | 2 +- src/TRDH_alg.jl | 2 +- src/TR_alg.jl | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl index afe70181..0df7c53f 100644 --- a/src/LMTR_alg.jl +++ b/src/LMTR_alg.jl @@ -508,6 +508,6 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(T), sqrt_ξ1_νInv) + set_residuals!(stats, zero(T), min(sqrt_ξ1_νInv, norm_s_cauchydν)) return stats end diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 93e0f6cd..65db1da4 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -499,6 +499,6 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(T), sqrt_ξ1_νInv) + set_residuals!(stats, zero(T), min(sqrt_ξ1_νInv, norm_s_cauchydν)) return stats end diff --git a/src/R2DH.jl b/src/R2DH.jl index 7310832c..2119dbe1 100644 --- a/src/R2DH.jl +++ b/src/R2DH.jl @@ -491,6 +491,6 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(eltype(xk)), sqrt_ξ_νInv) + set_residuals!(stats, zero(eltype(xk)), min(sqrt_ξ_νInv, norm_sdν)) return stats end diff --git a/src/R2N.jl b/src/R2N.jl index 96b607f2..22ba6e36 100644 --- a/src/R2N.jl +++ b/src/R2N.jl @@ -546,7 +546,7 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(eltype(xk)), sqrt_ξ1_νInv) + set_residuals!(stats, zero(eltype(xk)), min(sqrt_ξ1_νInv, norm_s_cauchydν)) return stats end diff --git a/src/R2_alg.jl b/src/R2_alg.jl index fce199a7..fc8d0b96 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -546,6 +546,6 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(eltype(xk)), sqrt_ξ_νInv) + set_residuals!(stats, zero(eltype(xk)), min(sqrt_ξ_νInv, norm_sdν)) return stats end diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl index 2bb0362c..1e165cf5 100644 --- a/src/TRDH_alg.jl +++ b/src/TRDH_alg.jl @@ -552,6 +552,6 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(eltype(xk)), sqrt_ξ_νInv) + set_residuals!(stats, zero(eltype(xk)), min(sqrt_ξ_νInv, norm_sdν)) return stats end diff --git a/src/TR_alg.jl b/src/TR_alg.jl index 8d7a3172..eeb9579b 100644 --- a/src/TR_alg.jl +++ b/src/TR_alg.jl @@ -558,5 +558,5 @@ function SolverCore.solve!( end set_solution!(stats, xk) - set_residuals!(stats, zero(eltype(xk)), sqrt_ξ1_νInv) + set_residuals!(stats, zero(eltype(xk)), min(sqrt_ξ1_νInv, norm_s_cauchydν)) end From 9453800dc4559a28d6eb926dbf0202ca2a430ec2 Mon Sep 17 00:00:00 2001 From: Maxence Gollier <134112149+MaxenceGollier@users.noreply.github.com> Date: Mon, 23 Feb 2026 09:57:59 -0500 Subject: [PATCH 7/9] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/LMTR_alg.jl | 6 +++--- src/LM_alg.jl | 2 +- src/R2DH.jl | 4 ++-- src/R2_alg.jl | 6 +++--- src/TRDH_alg.jl | 2 +- src/TR_alg.jl | 2 +- 6 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/LMTR_alg.jl b/src/LMTR_alg.jl index 0df7c53f..547e91d0 100644 --- a/src/LMTR_alg.jl +++ b/src/LMTR_alg.jl @@ -146,7 +146,7 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `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,)`. -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 or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ₁` is the Cauchy step. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -325,7 +325,7 @@ function SolverCore.solve!( atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) - + set_status!( stats, get_status( @@ -480,7 +480,7 @@ function SolverCore.solve!( norm_s_cauchy = norm(s) norm_s_cauchydν = norm_s_cauchy / ν - + solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) set_status!( diff --git a/src/LM_alg.jl b/src/LM_alg.jl index 65db1da4..186a3ada 100644 --- a/src/LM_alg.jl +++ b/src/LM_alg.jl @@ -146,7 +146,7 @@ For advanced usage, first define a solver "LMSolver" to preallocate the memory u - `subsolver = R2Solver`: the solver used to solve the subproblems. - `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,)`. -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 or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀₁‖/ν₀` where `sₖ₁` is the Cauchy step. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ₁‖/νₖ < atol_step + rtol_step*‖s₀₁‖/ν₀` where `sₖ₁` is the Cauchy step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. diff --git a/src/R2DH.jl b/src/R2DH.jl index 2119dbe1..f9e10ae7 100644 --- a/src/R2DH.jl +++ b/src/R2DH.jl @@ -360,7 +360,7 @@ function SolverCore.solve!( norm_s = norm(s) norm_sdν = norm_s / ν₁ atol_step += rtol_step * norm_sdν # make stopping test absolute and relative - + solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) set_status!( @@ -460,7 +460,7 @@ function SolverCore.solve!( sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν₁) : sqrt(-ξ / ν₁) (ξ < 0 && sqrt_ξ_νInv > neg_tol) && error("R2DH: prox-gradient step should produce a decrease but ξ = $(ξ)") - + norm_s = norm(s) norm_sdν = norm_s / ν₁ diff --git a/src/R2_alg.jl b/src/R2_alg.jl index fc8d0b96..516d97e1 100644 --- a/src/R2_alg.jl +++ b/src/R2_alg.jl @@ -160,7 +160,7 @@ For advanced usage, first define a solver "R2Solver" to preallocate the memory u - `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`; -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 or when `‖sₖ‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ` is the current step. +The algorithm stops either when `√(ξₖ/νₖ) < atol_decr + rtol_decr*√(ξ₀/ν₀) ` or `ξₖ < 0` and `√(-ξₖ/νₖ) < neg_tol` where ξₖ := f(xₖ) + h(xₖ) - φ(sₖ; xₖ) - ψ(sₖ; xₖ), and √(ξₖ/νₖ) is a stationarity measure or when `‖sₖ‖/νₖ < atol_step + rtol_step*‖s₀‖/ν₀` where `sₖ` is the current step. # Output The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. @@ -429,7 +429,7 @@ function SolverCore.solve!( atol_step += rtol_step * norm_sdν # make stopping test absolute and relative solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) - + set_status!( stats, get_status( @@ -520,7 +520,7 @@ function SolverCore.solve!( norm_sdν = norm_s / ν solved = (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) - + set_status!( stats, get_status( diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl index 1e165cf5..b612df96 100644 --- a/src/TRDH_alg.jl +++ b/src/TRDH_alg.jl @@ -411,7 +411,7 @@ function SolverCore.solve!( norm_sdν = norm_s / ν atol_step += rtol_step * norm_sdν # make stopping test absolute and relative - solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) | (norm_sdν ≤ atol_step) + solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) set_status!( stats, get_status( diff --git a/src/TR_alg.jl b/src/TR_alg.jl index eeb9579b..ec771916 100644 --- a/src/TR_alg.jl +++ b/src/TR_alg.jl @@ -363,7 +363,7 @@ function SolverCore.solve!( atol_step += rtol_step * norm_s_cauchydν # make stopping test absolute and relative solved = (ξ1 < 0 && sqrt_ξ1_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ1_νInv ≤ atol_decr) || (norm_s_cauchydν ≤ atol_step) - + set_status!( stats, get_status( From c75c8dc84d6b18f33a24825c6890cc59277f06a0 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Mon, 23 Feb 2026 10:08:33 -0500 Subject: [PATCH 8/9] solve variable scope --- src/TRDH_alg.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/TRDH_alg.jl b/src/TRDH_alg.jl index b612df96..612462f3 100644 --- a/src/TRDH_alg.jl +++ b/src/TRDH_alg.jl @@ -411,7 +411,8 @@ function SolverCore.solve!( norm_sdν = norm_s / ν atol_step += rtol_step * norm_sdν # make stopping test absolute and relative - solved = (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) + solved = reduce_TR ? (ξ1 < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ1 ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) : + (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv ≤ atol_decr) || (norm_sdν ≤ atol_step) set_status!( stats, get_status( From bf3de4133c05bb364f687ac1fa5b932452d39baf Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Mon, 23 Feb 2026 13:42:23 -0500 Subject: [PATCH 9/9] add tests --- test/runtests.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 2ceecf66..3d37177c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,12 +45,14 @@ for (mod, mod_name) ∈ ((x -> x, "exact"), (LSR1Model, "lsr1"), (LBFGSModel, "l @test length(out.solution) == bpdn.meta.nvar @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) out = solver(mod(bpdn), h, args..., options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @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 == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) end end end @@ -77,12 +79,14 @@ for (mod, mod_name) ∈ ((SpectralGradientModel, "spg"),) @test length(out.solution) == bpdn.meta.nvar @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) out = TRDH(mod(bpdn), h, χ, options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @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 == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) end end end @@ -107,12 +111,14 @@ for (mod, mod_name) ∈ ((LSR1Model, "lsr1"), (LBFGSModel, "lbfgs")) @test length(TR_out.solution) == bpdn.meta.nvar @test typeof(TR_out.dual_feas) == eltype(TR_out.solution) @test TR_out.status == :first_order + @test TR_out.step_status == (TR_out.iter > 0 ? :accepted : :unknown) TR_out = TR(mod(bpdn), h, NormL2(1.0), options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @test typeof(TR_out.solution) == typeof(bpdn.meta.x0) @test length(TR_out.solution) == bpdn.meta.nvar @test typeof(TR_out.dual_feas) == eltype(TR_out.solution) @test TR_out.status == :first_order + @test TR_out.step_status == (TR_out.iter > 0 ? :accepted : :unknown) end end end @@ -140,12 +146,14 @@ for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1"), (IndBallL0(10 * com @test length(out.solution) == bpdn.meta.nvar @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) out = solver(bpdn_nls, h, args..., options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @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 == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) end end end @@ -194,13 +202,14 @@ for (mod, mod_name) ∈ ( @test length(out.solution) == bpdn.meta.nvar @test typeof(out.dual_feas) == eltype(out.solution) @test out.status == :first_order + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) out = solver(mod(bpdn), h, options, x0 = x0, atol_decr = 0.0, rtol_decr = 0.0, atol_step = 1e-6, rtol_step = 1e-6) @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 == :first_order - + @test out.step_status == (out.iter > 0 ? :accepted : :unknown) end end end