diff --git a/src/R2DH.jl b/src/R2DH.jl new file mode 100644 index 00000000..1b48bc07 --- /dev/null +++ b/src/R2DH.jl @@ -0,0 +1,319 @@ +export R2DH + +""" + R2DH(nlp, h, options) + +A second-order quadratic regularization method for the problem + + min f(x) + h(x) + +where f: ℝⁿ → ℝ is C¹ and h: ℝⁿ → ℝ is lower semi-continuous, proper, and prox-bounded. + +About each iterate xₖ, a step sₖ is computed as a solution of + + min φ(s; xₖ) + ψ(s; xₖ) + +where φ(s ; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs + ½ sᵀ(Dₖ + σₖI)s is a quadratic approximation of f about xₖ, +ψ(s; xₖ) = h(xₖ + s), Dₖ is a diagonal Hessian approximation and σₖ > 0 is the regularization parameter. + +### Arguments + +* `nlp::AbstractDiagonalQNModel`: a smooth optimization problem +* `h`: a regularizer such as those defined in ProximalOperators +* `options::ROSolverOptions`: a structure containing algorithmic parameters + +### Keyword Arguments + +* `x0::AbstractVector`: an initial guess (in the first calling form: default = `nlp.meta.x0`) +* `selected::AbstractVector{<:Integer}`: subset of variables to which `h` is applied (default `1:length(x0)`). +* `D`: Diagonal quasi-Newton operator. +* `Mmonotone::Int`: number of previous values of the objective to consider for the non-monotone variant (default: 6). + +The objective and gradient of `nlp` will be accessed. + +### Return values +The value returned is a `GenericExecutionStats`, see `SolverCore.jl`. +""" +function R2DH( + nlp::AbstractDiagonalQNModel{R, S}, + h, + options::ROSolverOptions{R}; + kwargs..., + ) where {R <: Real, S} + kwargs_dict = Dict(kwargs...) + x0 = pop!(kwargs_dict, :x0, nlp.meta.x0) + xk, k, outdict = R2DH( + x -> obj(nlp, x), + (g, x) -> grad!(nlp, x, g), + h, + hess_op(nlp, x0), + options, + x0; + l_bound = nlp.meta.lvar, + u_bound = nlp.meta.uvar, + kwargs..., + ) + sqrt_ξ_νInv = outdict[:sqrt_ξ_νInv] + stats = GenericExecutionStats(nlp) + set_status!(stats, outdict[:status]) + set_solution!(stats, xk) + set_objective!(stats, outdict[:fk] + outdict[:hk]) + set_residuals!(stats, zero(eltype(xk)), sqrt_ξ_νInv) + set_iter!(stats, k) + set_time!(stats, outdict[:elapsed_time]) + set_solver_specific!(stats, :Fhist, outdict[:Fhist]) + set_solver_specific!(stats, :Hhist, outdict[:Hhist]) + set_solver_specific!(stats, :Time_hist, outdict[:Time_hist]) + set_solver_specific!(stats, :NonSmooth, outdict[:NonSmooth]) + set_solver_specific!(stats, :SubsolverCounter, outdict[:Chist]) + return stats + end + +""" + R2DH(f, ∇f!, h, options, x0) + +A second calling form for `R2DH` where the objective and gradient are passed as arguments. + +### Arguments + +* `f::Function`: the objective function +* `∇f!::Function`: the gradient function +* `h`: a regularizer such as those defined in ProximalOperators +* `D`: Diagonal quasi-Newton operator. +* `options::ROSolverOptions`: a structure containing algorithmic parameters +* `x0::AbstractVector`: an initial guess + +### Keyword Arguments + +* `Mmonotone::Int`: number of previous values of the objective to consider for the non-monotone variant (default: 6). +* `selected::AbstractVector{<:Integer}`: subset of variables to which `h` is applied (default `1:length(x0)`). + +### Return values + +* `xk`: the final iterate +* `k`: the number of iterations +* `outdict`: a dictionary containing the following fields: + * `Fhist`: an array with the history of values of the smooth objective + * `Hhist`: an array with the history of values of the nonsmooth objective + * `Time_hist`: an array with the history of elapsed times + * `Chist`: an array with the history of number of inner iterations + * `NonSmooth`: the nonsmooth term + * `status`: the status of the solver either `:first_order`, `:max_iter`, `:max_time` or `:exception` + * `fk`: the value of the smooth objective at the final iterate + * `hk`: the value of the nonsmooth objective at the final iterate + * `sqrt_ξ_νInv`: the square root of the ratio of the nonsmooth term to the regularization parameter + * `elapsed_time`: the elapsed time +""" +function R2DH( + f::F, + ∇f!::G, + h::H, + D::DQN, + options::ROSolverOptions{R}, + x0::AbstractVector{R}; + Mmonotone::Int = 6, + selected::AbstractVector{<:Integer} = 1:length(x0), + kwargs..., +) where {F <: Function, G <: Function, H, R <: Real, DQN <: AbstractDiagonalQuasiNewtonOperator} + start_time = time() + elapsed_time = 0.0 + ϵ = options.ϵa + ϵr = options.ϵr + neg_tol = options.neg_tol + verbose = options.verbose + maxIter = options.maxIter + maxTime = options.maxTime + σmin = options.σmin + σk = options.σk + η1 = options.η1 + η2 = options.η2 + ν = options.ν + γ = options.γ + θ = options.θ + + local l_bound, u_bound + has_bnds = false + for (key, val) in kwargs + if key == :l_bound + l_bound = val + has_bnds = has_bnds || any(l_bound .!= R(-Inf)) + elseif key == :u_bound + u_bound = val + has_bnds = has_bnds || any(u_bound .!= R(Inf)) + end + end + + if verbose == 0 + ptf = Inf + elseif verbose == 1 + ptf = round(maxIter / 10) + elseif verbose == 2 + ptf = round(maxIter / 100) + else + ptf = 1 + end + + # initialize parameters + xk = copy(x0) + hk = h(xk[selected]) + if hk == Inf + verbose > 0 && @info "R2DH: finding initial guess where nonsmooth term is finite" + prox!(xk, h, x0, one(eltype(x0))) + hk = h(xk[selected]) + hk < Inf || error("prox computation must be erroneous") + verbose > 0 && @debug "R2DH: found point where h has value" hk + end + hk == -Inf && error("nonsmooth term is not proper") + + xkn = similar(xk) + s = zero(xk) + ψ = has_bnds ? shifted(h, xk, l_bound - xk, u_bound - xk, selected) : shifted(h, xk) + + Fobj_hist = zeros(maxIter+1) + Hobj_hist = zeros(maxIter+1) + time_hist = zeros(maxIter+1) + FHobj_hist = fill!(Vector{R}(undef, Mmonotone - 1), R(-Inf)) + Complex_hist = zeros(Int, maxIter+1) + if verbose > 0 + #! format: off + @info @sprintf "%6s %8s %8s %7s %8s %7s %7s %7s %1s" "iter" "f(x)" "h(x)" "√(ξ/ν)" "ρ" "σ" "‖x‖" "‖s‖" "" + #! format: off + end + + local ξ + k = 0 + + fk = f(xk) + ∇fk = similar(xk) + ∇f!(∇fk, xk) + ∇fk⁻ = copy(∇fk) + spectral_test = isa(D, SpectralGradient) + Dkσk = D.d .+ σk + DNorm = norm(D.d, Inf) + + ν = 1 / ((DNorm + σk) * (1 + θ)) + mν∇fk = -ν * ∇fk + sqrt_ξ_νInv = one(R) + + optimal = false + tired = maxIter > 0 && k ≥ maxIter || elapsed_time > maxTime + + while !(optimal || tired) + # model with diagonal hessian + φ(d) = ∇fk' * d + (d' * (Dkσk .* d)) / 2 + mk(d) = φ(d) + ψ(d) + + if spectral_test + prox!(s, ψ, mν∇fk, ν) + else + iprox!(s, ψ, ∇fk, Dkσk) + end + mks = mk(s) + + k = k + 1 + elapsed_time = time() - start_time + Fobj_hist[k] = fk + Hobj_hist[k] = hk + time_hist[k] = elapsed_time + Mmonotone > 1 && (FHobj_hist[mod(k-1, Mmonotone - 1) + 1] = fk + hk) + + Complex_hist[k] += 1 + xkn .= xk .+ s + fkn = f(xkn) + hkn = h(xkn[selected]) + hkn == -Inf && error("nonsmooth term is not proper") + + fhmax = Mmonotone > 1 ? maximum(FHobj_hist) : fk + hk + Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps() + Δmod = fhmax - (fk + mks) + max(1, abs(hk)) * 10 * eps() + ξ = hk - mks + max(1, abs(hk)) * 10 * eps() + sqrt_ξ_νInv = ξ ≥ 0 ? sqrt(ξ / ν) : sqrt(-ξ / ν) + + if ξ ≥ 0 && k == 1 + ϵ += ϵr * sqrt_ξ_νInv # make stopping test absolute and relative + end + + if (ξ < 0 && sqrt_ξ_νInv ≤ neg_tol) || (ξ ≥ 0 && sqrt_ξ_νInv < ϵ) + # the current xk is approximately first-order stationary + optimal = true + continue + end + + if (ξ ≤ 0 || isnan(ξ)) + error("R2DH: failed to compute a step: ξ = $ξ") + end + + ρk = Δobj / Δmod + + σ_stat = (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=") + + if (verbose > 0) && ((k % ptf == 0) || (k == 1)) + #! format: off + @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e %7.1e %7.1e %1s" k fk hk sqrt_ξ_νInv ρk σk norm(xk) norm(s) σ_stat + #! format: on + end + + if η2 ≤ ρk < Inf + σk = max(σk / γ, σmin) + end + + if η1 ≤ ρk < Inf + xk .= xkn + has_bnds && set_bounds!(ψ, l_bound - xk, u_bound - xk) + fk = fkn + hk = hkn + shift!(ψ, xk) + ∇f!(∇fk, xk) + push!(D, s, ∇fk - ∇fk⁻) # update QN operator + DNorm = norm(D.d, Inf) + ∇fk⁻ .= ∇fk + end + + if ρk < η1 || ρk == Inf + σk = σk * γ + end + + Dkσk .= D.d .+ σk + ν = 1 / ((DNorm + σk) * (1 + θ)) + + tired = maxIter > 0 && k ≥ maxIter + if !tired + @. mν∇fk = -ν * ∇fk + end + end + + if verbose > 0 + if k == 1 + @info @sprintf "%6d %8.1e %8.1e" k fk hk + elseif optimal + #! format: off + @info @sprintf "%6d %8.1e %8.1e %7.1e %8s %7.1e %7.1e %7.1e" k fk hk sqrt_ξ_νInv "" σk norm(xk) norm(s) + #! format: on + @info "R2DH: terminating with √(ξ/ν) = $(sqrt_ξ_νInv))" + end + end + + status = if optimal + :first_order + elseif elapsed_time > maxTime + :max_time + elseif tired + :max_iter + else + :exception + end + outdict = Dict( + :Fhist => Fobj_hist[1:k], + :Hhist => Hobj_hist[1:k], + :Time_hist => time_hist[1:k], + :Chist => Complex_hist[1:k], + :NonSmooth => h, + :status => status, + :fk => fk, + :hk => hk, + :sqrt_ξ_νInv => sqrt_ξ_νInv, + :elapsed_time => elapsed_time, + ) + + return xk, k, outdict +end diff --git a/src/R2N.jl b/src/R2N.jl new file mode 100644 index 00000000..f914372a --- /dev/null +++ b/src/R2N.jl @@ -0,0 +1,294 @@ +export R2N + +""" +R2N(nlp, h, χ, options; kwargs...) + +A regularized quasi-Newton method for the problem + + min f(x) + h(x) + +where f: ℝⁿ → ℝ is C¹ and h: ℝⁿ → ℝ is lower semi-continuous, proper and prox-bounded. + +About each iterate xₖ, a step sₖ is computed as an approximate solution of + + min φ(s; xₖ) + ½ σₖ ‖s‖² + ψ(s; xₖ) + +where φ(s; xₖ) = f(xₖ) + ∇f(xₖ)ᵀs + ½ sᵀBₖs is a quadratic approximation of f about xₖ, +ψ(s; xₖ) = h(xₖ + s) and σₖ > 0 is the regularization parameter. +The subproblem is solved inexactly by way of a first-order method such as the proximal-gradient +method or the quadratic regularization method. + +### Arguments + +* `nlp::AbstractNLPModel`: a smooth optimization problem +* `h`: a regularizer such as those defined in ProximalOperators +* `options::ROSolverOptions`: a structure containing algorithmic parameters. + +The objective, gradient and Hessian of `nlp` will be accessed. +The Hessian is accessed as an abstract operator and need not be the exact Hessian. + +### Keyword arguments + +* `x0::AbstractVector`: an initial guess (default: `nlp.meta.x0`) +* `subsolver_logger::AbstractLogger`: a logger to pass to the subproblem solver (default: the null logger) +* `subsolver`: the procedure used to compute a step (`R2DH`, `R2` or `PG`) +* `subsolver_options::ROSolverOptions`: default options to pass to the subsolver (default: all defaut options) +* `Mmonotone::Int`: number of previous values of the objective to consider for the non-monotone variant (default: 1). +* `selected::AbstractVector{<:Integer}`: subset of variables to which `h` is applied (default `1:nlp.meta.nvar`). + +### Return values + +* `xk`: the final iterate +* `Fobj_hist`: an array with the history of values of the smooth objective +* `Hobj_hist`: an array with the history of values of the nonsmooth objective +* `Complex_hist`: an array with the history of number of inner iterations. +""" +function R2N( + f::AbstractNLPModel, + h::H, + options::ROSolverOptions{R}; + x0::AbstractVector = f.meta.x0, + subsolver_logger::Logging.AbstractLogger = Logging.NullLogger(), + subsolver = R2, + subsolver_options = ROSolverOptions(ϵa = options.ϵa), + Mmonotone::Int = 1, + selected::AbstractVector{<:Integer} = 1:(f.meta.nvar), +) where {H, R} + start_time = time() + elapsed_time = 0.0 + # initialize passed options + ϵ = options.ϵa + ϵ_subsolver_init = subsolver_options.ϵa + ϵ_subsolver = copy(ϵ_subsolver_init) + ϵr = options.ϵr + Δk = options.Δk + verbose = options.verbose + maxIter = options.maxIter + maxTime = options.maxTime + η1 = options.η1 + η2 = options.η2 + γ = options.γ + θ = options.θ + σmin = options.σmin + α = options.α + β = options.β + σk = options.σk + + # store initial values of the subsolver_options fields that will be modified + ν_subsolver = subsolver_options.ν + ϵa_subsolver = subsolver_options.ϵa + + local l_bound, u_bound + if has_bounds(f) + l_bound = f.meta.lvar + u_bound = f.meta.uvar + end + + if verbose == 0 + ptf = Inf + elseif verbose == 1 + ptf = round(maxIter / 10) + elseif verbose == 2 + ptf = round(maxIter / 100) + else + ptf = 1 + end + + # initialize parameters + xk = copy(x0) + hk = h(xk[selected]) + if hk == Inf + verbose > 0 && @info "R2N: finding initial guess where nonsmooth term is finite" + prox!(xk, h, x0, one(eltype(x0))) + hk = h(xk[selected]) + hk < Inf || error("prox computation must be erroneous") + verbose > 0 && @debug "R2N: found point where h has value" hk + end + hk == -Inf && error("nonsmooth term is not proper") + + xkn = similar(xk) + s = zero(xk) + ψ = has_bounds(f) ? shifted(h, xk, l_bound - xk, u_bound - xk, selected) : shifted(h, xk) + + Fobj_hist = zeros(maxIter) + Hobj_hist = zeros(maxIter) + FHobj_hist = fill!(Vector{R}(undef, Mmonotone - 1), R(-Inf)) + Complex_hist = zeros(Int, maxIter) + 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‖" "‖Bₖ‖" "R2N" + #! format: on + end + + # main algorithm initialization + + local ξ1 + k = 0 + + fk = obj(f, xk) + ∇fk = grad(f, xk) + ∇fk⁻ = copy(∇fk) + + quasiNewtTest = isa(f, QuasiNewtonModel) + Bk = hess_op(f, xk) + + λmax = opnorm(Bk) + νInv = (1 + θ) * (σk + λmax) + sqrt_ξ1_νInv = one(R) + + optimal = false + tired = k ≥ maxIter || elapsed_time > maxTime + + while !(optimal || tired) + k = k + 1 + elapsed_time = time() - start_time + Fobj_hist[k] = fk + Hobj_hist[k] = hk + Mmonotone > 1 && (FHobj_hist[mod(k-1, Mmonotone - 1) + 1] = fk + hk) + + # model for first prox-gradient step and ξ1 + φ1(d) = ∇fk' * d + mk1(d) = φ1(d) + ψ(d) + + # model for subsequent prox-gradient steps and ξ + φ(d) = ∇fk' * d + dot(d, Bk * d) / 2 + σk * dot(d, d) / 2 + + ∇φ!(g, d) = begin + mul!(g, Bk, d) + g .+= ∇fk + g .+= σk * d + g + end + + mk(d) = φ(d) + ψ(d) + + # take first proximal gradient step s1 and see if current xk is nearly stationary + # s1 minimizes φ1(s) + ‖s‖² / 2 / ν + ψ(s) ⟺ s1 ∈ prox{νψ}(-ν∇φ1(0)). + + subsolver_options.ν = 1 / νInv + prox!(s, ψ, -subsolver_options.ν * ∇fk, subsolver_options.ν) + ξ1 = hk - mk1(s) + max(1, abs(hk)) * 10 * eps() + ξ1 > 0 || error("R2N: 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_νInv + ϵ += ϵ_increment # make stopping test absolute and relative + ϵ_subsolver += ϵ_increment + end + + if sqrt_ξ1_νInv < ϵ + # the current xk is approximately first-order stationary + optimal = true + continue + end + s1 = copy(s) + + subsolver_options.ϵa = k == 1 ? 1.0e-3 : min(sqrt_ξ1_νInv ^ (1.5) , sqrt_ξ1_νInv * 1e-3) + verbose > 0 && @debug "setting inner stopping tolerance to" subsolver_options.optTol + subsolver_options.σk = σk + subsolver_args = subsolver == R2DH ? (SpectralGradient(νInv, f.meta.nvar),) : () + s, iter, _ = with_logger(subsolver_logger) do + subsolver(φ, ∇φ!, ψ, subsolver_args..., subsolver_options, s) + end + + if norm(s) > β * norm(s1) + s .= s1 + end + + # restore initial subsolver_options.ϵa here so that subsolver_options.ϵa + # is not modified if there is an error + subsolver_options.ν = ν_subsolver + subsolver_options.ϵa = ϵ_subsolver_init + subsolver_options.σk = σk + Complex_hist[k] = iter + + xkn .= xk .+ s + fkn = obj(f, xkn) + hkn = h(xkn[selected]) + hkn == -Inf && error("nonsmooth term is not proper") + mks = mk(s) + + fhmax = Mmonotone > 1 ? maximum(FHobj_hist) : fk + hk + Δobj = fhmax - (fkn + hkn) + max(1, abs(fhmax)) * 10 * eps() + Δmod = fhmax - (fk + mks) + max(1, abs(fhmax)) * 10 * eps() + ξ = hk - mks + max(1, abs(hk)) * 10 * eps() + + if (ξ ≤ 0 || isnan(ξ)) + error("R2N: failed to compute a step: ξ = $ξ") + end + + ρk = Δobj / Δmod + + R2N_stat = (η2 ≤ ρk < Inf) ? "↘" : (ρk < η1 ? "↗" : "=") + + if (verbose > 0) && ((k % ptf == 0) || (k == 1)) + #! 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_νInv sqrt(ξ1) ρk σk norm(xk) norm(s) λmax R2N_stat + #! format: off + end + + if η2 ≤ ρk < Inf + σk = max(σk/γ, σmin) + end + + if η1 ≤ ρk < Inf + xk .= xkn + has_bounds(f) && set_bounds!(ψ, l_bound - xk, u_bound - xk) + + #update functions + fk = fkn + hk = hkn + + # update gradient & Hessian + shift!(ψ, xk) + ∇fk = grad(f, xk) + if quasiNewtTest + push!(f, s, ∇fk - ∇fk⁻) + end + Bk = hess_op(f, xk) + λmax = opnorm(Bk) + ∇fk⁻ .= ∇fk + end + + if ρk < η1 || ρk == Inf + σk = σk * γ + end + νInv = (1 + θ) * (σk + λmax) + tired = k ≥ maxIter || elapsed_time > maxTime + end + + if verbose > 0 + if k == 1 + @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_νInv sqrt(ξ1) "" σk norm(xk) norm(s) λmax + #! format: on + @info "R2N: terminating with √(ξ1/ν) = $(sqrt_ξ1_νInv)" + end + end + + status = if optimal + :first_order + elseif elapsed_time > maxTime + :max_time + elseif tired + :max_iter + else + :exception + end + + stats = GenericExecutionStats(f) + set_status!(stats, status) + set_solution!(stats, xk) + set_objective!(stats, fk + hk) + set_residuals!(stats, zero(eltype(xk)), sqrt_ξ1_νInv) + set_iter!(stats, k) + set_time!(stats, elapsed_time) + set_solver_specific!(stats, :Fhist, Fobj_hist[1:k]) + set_solver_specific!(stats, :Hhist, Hobj_hist[1:k]) + set_solver_specific!(stats, :NonSmooth, h) + set_solver_specific!(stats, :SubsolverCounter, Complex_hist[1:k]) + return stats +end \ No newline at end of file diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl index 22c1c814..5983cab9 100644 --- a/src/RegularizedOptimization.jl +++ b/src/RegularizedOptimization.jl @@ -20,5 +20,7 @@ include("TRDH_alg.jl") include("R2_alg.jl") include("LM_alg.jl") include("LMTR_alg.jl") +include("R2DH.jl") +include("R2N.jl") end # module RegularizedOptimization diff --git a/src/input_struct.jl b/src/input_struct.jl index 2e9937ce..7242b90a 100644 --- a/src/input_struct.jl +++ b/src/input_struct.jl @@ -9,6 +9,7 @@ mutable struct ROSolverOptions{R} maxIter::Int # maximum amount of inner iterations maxTime::Float64 #maximum time allotted to the algorithm in s σmin::R # minimum σk allowed for LM/R2 method + σk::R # initial σk η1::R # step acceptance threshold η2::R # trust-region increase threshold α::R # νk Δ^{-1} parameter @@ -27,6 +28,7 @@ mutable struct ROSolverOptions{R} maxIter::Int = 10000, maxTime::Float64 = 3600.0, σmin::R = eps(R), + σk::R = eps(R)^(1 / 5), η1::R = √√eps(R), η2::R = R(0.9), α::R = 1 / eps(R), @@ -44,6 +46,7 @@ mutable struct ROSolverOptions{R} @assert maxIter ≥ 0 @assert maxTime ≥ 0 @assert σmin ≥ 0 + @assert σk ≥ 0 @assert 0 < η1 < η2 < 1 @assert α > 0 @assert ν > 0 @@ -59,6 +62,7 @@ mutable struct ROSolverOptions{R} maxIter, maxTime, σmin, + σk, η1, η2, α, diff --git a/test/runtests.jl b/test/runtests.jl index 39f3a9a0..b496c29e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -134,5 +134,33 @@ for (h, h_name) ∈ ((NormL1(λ), "l1"),) end end +R2N_R2DH(args...; kwargs...) = R2N(args...; subsolver = R2DH, kwargs...) +for (mod, mod_name) ∈ ((SpectralGradientModel, "spg"), (DiagonalPSBModel, "psb"), (LSR1Model, "lsr1"), (LBFGSModel, "lbfgs")) + for (h, h_name) ∈ ((NormL0(λ), "l0"), (NormL1(λ), "l1")) + for solver_sym ∈ (:R2DH, :R2N, :R2N_R2DH) + solver_sym ∈ (:R2N,:R2N_R2DH) && mod_name ∈ ("spg", "psb") && continue + solver_sym == :R2DH && mod_name != "spg" && continue + solver_sym == :R2N_R2DH && h_name == "l1" && continue # this test seems to fail because s seems to be equal to zeros within the subsolver + solver_name = string(solver_sym) + solver = eval(solver_sym) + @testset "bpdn-$(mod_name)-$(solver_name)-$(h_name)" begin + x0 = zeros(bpdn.meta.nvar) + out = solver(mod(bpdn), h, options, x0 = x0) + @test typeof(out.solution) == typeof(bpdn.meta.x0) + @test length(out.solution) == bpdn.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, out.solution) == out.solver_specific[:Fhist][end] + @test h(out.solution) == out.solver_specific[:Hhist][end] + @test out.status == :first_order + end + end + end +end + include("test_bounds.jl") include("test_allocs.jl")