From a8392438ddbcb457ce9b67f2868688af7645d021 Mon Sep 17 00:00:00 2001 From: Maxence Gollier Date: Sat, 6 Sep 2025 13:43:04 -0400 Subject: [PATCH] remove dead code --- src/Fista_alg.jl | 232 ------------------------- src/PG_alg.jl | 307 --------------------------------- src/RegularizedOptimization.jl | 3 - src/splitting.jl | 213 ----------------------- 4 files changed, 755 deletions(-) delete mode 100644 src/Fista_alg.jl delete mode 100644 src/PG_alg.jl delete mode 100644 src/splitting.jl diff --git a/src/Fista_alg.jl b/src/Fista_alg.jl deleted file mode 100644 index 33bacd72..00000000 --- a/src/Fista_alg.jl +++ /dev/null @@ -1,232 +0,0 @@ -export FISTA, FISTAD - -""" - FISTA for - min_x ϕ(x) = f(x) + g(x), with f(x) cvx and β-smooth, g(x) closed cvx - - Input: - f: function handle that returns f(x) and ∇f(x) - h: function handle that returns g(x) - s: initial point - proxG: function handle that calculates prox_{νg} - options: see descentopts.jl - Output: - s⁺: s update - s : s^(k-1) - his : function history - feval : number of function evals (total objective) -""" -function FISTA(nlp::AbstractNLPModel, args...; kwargs...) - kwargs_dict = Dict(kwargs...) - x0 = pop!(kwargs_dict, :x0, nlp.meta.x0) - xk, k, outdict = FISTA(x -> obj(nlp, x), (g, x) -> grad!(nlp, x, g), args..., x0; kwargs_dict...) - - return GenericExecutionStats( - outdict[:status], - nlp, - solution = xk, - objective = outdict[:fk] + outdict[:hk], - dual_feas = sqrt(outdict[:ξ]), - iter = k, - elapsed_time = outdict[:elapsed_time], - solver_specific = Dict( - :Fhist => outdict[:Fhist], - :Hhist => outdict[:Hhist], - :NonSmooth => outdict[:NonSmooth], - :SubsolverCounter => outdict[:Chist], - ), - ) -end - -function FISTA( - f::F, - ∇f!::G, - h::H, - options::ROSolverOptions, - x0::AbstractVector, -) where {F <: Function, G <: Function, H} - start_time = time() - elapsed_time = 0.0 - ϵ = options.ϵ - maxIter = options.maxIter - maxTime = options.maxTime - ν = options.ν - verbose = options.verbose - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - - #Problem Initialize - xk = copy(x0) - y = similar(xk) - ∇fk = zero(xk) - ∇fkn = similar(∇fk) - xkn = zero(xk) - fstep = xk .- ν .* ∇fk - Fobj_hist = zeros(maxIter) - Hobj_hist = zeros(maxIter) - Complex_hist = zeros(Int, maxIter) - - #initialize parameters - t = 1.0 - local ξ - k = 0 - ∇f!(∇fk, xk) #objInner/ quadratic model - fk = f(xk) - hk = h(xk) - - optimal = false - tired = k ≥ maxIter || elapsed_time ≥ maxTime - - if options.verbose != 0 - @info @sprintf "%6s %8s %8s %7s %8s %7s" "iter" "f(x)" "h(x)" "‖∂ϕ‖" "ν" "‖x‖" - end - - while !(optimal || tired) - k = k + 1 - elapsed_time = time() - start_time - Fobj_hist[k] = fk - Hobj_hist[k] = hk - Complex_hist[k] += 1 - - ∇fkn .= ∇fk - xkn .= xk - fstep .= xk .- ν .* ∇fk - prox!(xk, h, fstep, ν) - - #update step - t⁻ = t - t = 0.5 * (1.0 + sqrt(1.0 + 4.0 * t⁻^2)) - - #update y - y .= xk .+ ((t⁻ - 1.0) / t) .* (xk .- xkn) - - ∇f!(∇fk, xk) - fk = f(xk) - hk = h(xk) - - k += 1 - ∇fkn .= ∇fk .- ∇fkn .- (xk .- xkn) ./ ν - ξ = norm(∇fkn) - optimal = ξ < ϵ - tired = k ≥ maxIter || elapsed_time > maxTime - - if (verbose > 0) && (k % ptf == 0) - @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e " k fk hk ξ ν norm(xk) - 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], - :Chist => Complex_hist[1:k], - :NonSmooth => h, - :status => status, - :fk => fk, - :hk => hk, - :ξ => ξ, - :elapsed_time => elapsed_time, - ) - - return xk, k, outdict -end - -#enforces strict descent for FISTA -function FISTAD(f, ∇f, h, options; x0::AbstractVector = f.meta.x0) - ϵ = options.ϵ - maxIter = options.maxIter - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - - #Problem Initialize - ν = options.ν - v = deepcopy(x) - x⁺ = zero(x) - Fobj_hist = zeros(maxIter) - Hobj_hist = zeros(maxIter) - Complex_hist = zeros(Int, maxIter) - - #initialize parameters - t = 1.0 - # Iteration set up - k = 1 - - #do iterations - y = (1.0 - t) * x + t * v - g = ∇f(y) - fk = f(y) - hk = h(y) - - optimal = false - tired = k ≥ maxIter - - if options.verbose != 0 - @info @sprintf "%6s %8s %8s %7s %8s %7s" "iter" "f(x)" "h(x)" "‖∂ϕ‖" "ν" "‖x‖" - end - - while !(optimal || tired) - copy!(x, x⁺) - gold = g - - #complete prox step - u = ShiftedProximalOperators.prox(h, y - ν * g, ν) - - if f(u) ≤ f #this does not work - x⁺ = u - else - x⁺ = x - end - - #update step - # t⁻ = t - # t = R(0.5)*(R(1.0) + sqrt(R(1.0)+R(4.0)*t⁻^2)) - t = 2 / (k + 1) - - #update y - # v = s⁺ + ((t⁻ - R(1.0))/t)*(s⁺-s) - v = x⁺ + (1.0 / t) * (u - s⁺) - y = (1.0 - t) * x⁺ + t * v #I think this shold be s⁺ since it's at the end of the loop - - #update parameters - g = ∇f(y) - f = f(y) - - #check convergence - err = norm(g - gold - (x⁺ - x) / ν) - k += 1 - optimal = err < ϵ - tired = k ≥ maxIter - - k % ptf == 0 && @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e " k fk hk err ν norm(xk) - end - return x⁺, - k, - Fobj_hist[Fobj_hist .!= 0], - Hobj_hist[Fobj_hist .!= 0], - Complex_hist[Complex_hist .!= 0] -end diff --git a/src/PG_alg.jl b/src/PG_alg.jl deleted file mode 100644 index 6c66212e..00000000 --- a/src/PG_alg.jl +++ /dev/null @@ -1,307 +0,0 @@ -export PG, PGLnsch, PGΔ, PGE - -""" -Proximal Gradient Descent for - - min_x ϕ(x) = f(x) + g(x), with f(x) β-smooth, g(x) closed, lsc - -Input: - f: function handle that returns f(x) and ∇f(x) - h: function handle that returns g(x) - s: initial point - proxG: function handle that calculates prox_{νg} - options: see descentopts.jl -Output: - s⁺: s update - s : s^(k-1) - his : function history - feval : number of function evals (total objective ) -""" -function PG(nlp::AbstractNLPModel, args...; kwargs...) - kwargs_dict = Dict(kwargs...) - x0 = pop!(kwargs_dict, :x0, nlp.meta.x0) - xk, k, outdict = PG(x -> obj(nlp, x), (g, x) -> grad!(nlp, x, g), args..., x0; kwargs_dict...) - - return GenericExecutionStats( - outdict[:status], - nlp, - solution = xk, - objective = outdict[:fk] + outdict[:hk], - dual_feas = sqrt(outdict[:ξ]), - iter = k, - elapsed_time = outdict[:elapsed_time], - solver_specific = Dict( - :Fhist => outdict[:Fhist], - :Hhist => outdict[:Hhist], - :NonSmooth => outdict[:NonSmooth], - :SubsolverCounter => outdict[:Chist], - ), - ) -end - -function PG( - f::F, - ∇f!::G, - h::H, - options::ROSolverOptions, - x0::AbstractVector, -) where {F <: Function, G <: Function, H} - start_time = time() - elapsed_time = 0.0 - ϵ = options.ϵa - maxIter = options.maxIter - maxTime = options.maxTime - ν = options.ν - verbose = options.verbose - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - - #Problem Initialize - xk = copy(x0) - Fobj_hist = zeros(maxIter) - Hobj_hist = zeros(maxIter) - Complex_hist = zeros(Int, maxIter) - - # Iteration set up - ∇fk = similar(xk) - ∇f!(∇fk, xk) - fk = f(xk) - hk = h(xk) - ∇fkn = similar(∇fk) - xkn = similar(x0) - fstep = xk .- ν .* ∇fk - - #do iterations - local ξ - k = 0 - optimal = false - tired = k ≥ maxIter || elapsed_time ≥ maxTime - - if options.verbose != 0 - @info @sprintf "%6s %8s %8s %7s %8s %7s" "iter" "f(x)" "h(x)" "‖∂ϕ‖" "ν" "‖x‖" - end - - while !(optimal || tired) - elapsed_time = time() - start_time - Fobj_hist[k] = fk - Hobj_hist[k] = hk - Complex_hist[k] += 1 - - ∇fkn .= ∇fk - xkn .= xk - fstep .= xk .- ν .* ∇fk - prox!(xk, h, fstep, ν) - - ∇f!(∇fk, xk) - fk = f(xk) - hk = h(xk) - - k += 1 - ∇fkn .= ∇fk .- ∇fkn .- (xk .- xkn) ./ ν - ξ = norm(∇fkn) - optimal = ξ < ϵ - tired = k ≥ maxIter || elapsed_time > maxTime - - if (verbose > 0) && (k % ptf == 0) - @info @sprintf "%6d %8.1e %8.1e %7.1e %8.1e %7.1e " k fk hk ξ ν norm(xk) - 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], - :Chist => Complex_hist[1:k], - :NonSmooth => h, - :status => status, - :fk => fk, - :hk => hk, - :ξ => ξ, - :elapsed_time => elapsed_time, - ) - - return xk, k, outdict -end - -function PGΔ(f, ∇f, h, options; x::AbstractVector = f.meta.x0) - ε = options.optTol - maxIter = options.maxIter - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - #Problem Initialize - ν = options.ν - p = options.p - fDec = options.fDec - - k = 1 - feval = 1 - x⁺ = deepcopy(x) - - # Iteration set up - g = ∇f(x⁺) #objInner/ quadratic model - fk = f(x⁺) - - #do iterations - optimal = false - FD = false - tired = k ≥ maxIter - - while !(optimal || tired || FD) - gold = g - fold = f - x = x⁺ - - x⁺ = ShiftedProximalOperators.prox(h, x - ν * g, ν) - # update function info - g = ∇f(x⁺) - f = f(x⁺) - - feval += 1 - k += 1 - err = norm(g - gold - (x⁺ - x) / ν) - optimal = err < ε - tired = k ≥ maxIter - - k % ptf == 0 && @info @sprintf "%4d ‖xᵏ⁺¹ - xᵏ‖=%1.5e ν = %1.5e" k err ν - - Difff = fold + h(x) - f - h(x⁺) # these do not work - FD = abs(Difff) < p * norm(fDec) - end - return x⁺, feval -end - -function PGE(f, h, s, options) - ε = options.optTol - maxIter = options.maxIter - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - #Problem Initialize - ν = options.ν - p = options.p - fDec = options.fDec - - k = 1 - feval = 1 - s⁺ = deepcopy(s) - - # Iteration set up - g = ∇f(s⁺) #objInner/ quadratic model - - #do iterations - FD = false - optimal = false - tired = k ≥ maxIter - - #do iterations - while !(optimal || tired || FD) - gold = g - s = s⁺ - - ν = min(g' * g / (g' * Bk * g), ν) #no BK, will not work - #prox step - s⁺ = ShiftedProximalOperators.prox(h, s - ν * g, ν) - # update function info - g = ∇f(s⁺) - f = f(s⁺) - - feval += 1 - k += 1 - err = norm(g - gold - (s⁺ - s) / ν) - optimal = err < ε - tired = k ≥ maxIter - - k % ptf == 0 && @info @sprintf "%4d ‖xᵏ⁺¹ - xᵏ‖=%1.5e ν = %1.5e" k err ν - - Difff = fold + h(s) - f - h(s⁺) # these do not work - FD = abs(Difff) < p * norm(fDec) - end - return s⁺, feval -end - -function PGLnsch(f, ∇f, h, s, options) - ε = options.optTol - maxIter = options.maxIter - - if options.verbose == 0 - ptf = Inf - elseif options.verbose == 1 - ptf = round(maxIter / 10) - elseif options.verbose == 2 - ptf = round(maxIter / 100) - else - ptf = 1 - end - - #Problem Initialize - p = options.p - ν₀ = options.ν - k = 1 - s⁺ = deepcopy(s) - - # Iteration set up - feval = 1 - g = ∇f(s⁺) #objInner/ quadratic model - fk = f(s⁺) - - #do iterations - optimal = false - tired = k ≥ maxIter - - while !(optimal || tired) - gold = g - s = s⁺ - - s⁺ = ShiftedProximalOperators.prox(h, s - ν * g, ν) - #linesearch but we don't pass in f? - while f(s⁺) ≥ fk + g' * (s⁺ - s) + 1 / (ν * 2) * norm(s⁺ - s)^2 - ν *= p * ν - s⁺ = prox(h, s - ν * g, ν) - feval += 1 - end - # update function info - g = ∇f(s⁺) #objInner/ quadratic model - fk = f(s⁺) - - feval += 1 - k += 1 - err = norm(g - gold - (s⁺ - s) / ν) - optimal = err < ε - tired = k ≥ maxIter - - k % ptf == 0 && @info @sprintf "%4d ‖xᵏ⁺¹ - xᵏ‖=%1.5e ν = %1.5e" k err ν - ν = ν₀ - end - return s⁺, feval -end diff --git a/src/RegularizedOptimization.jl b/src/RegularizedOptimization.jl index 01778e86..fe13c4b9 100644 --- a/src/RegularizedOptimization.jl +++ b/src/RegularizedOptimization.jl @@ -13,9 +13,6 @@ using Percival: AugLagModel, update_y!, update_μ! include("utils.jl") include("input_struct.jl") -include("PG_alg.jl") -include("Fista_alg.jl") -include("splitting.jl") include("TR_alg.jl") include("TRDH_alg.jl") include("R2_alg.jl") diff --git a/src/splitting.jl b/src/splitting.jl deleted file mode 100644 index acf51a92..00000000 --- a/src/splitting.jl +++ /dev/null @@ -1,213 +0,0 @@ -export prox_split_1w, prox_split_2w - -"""Solves descent direction s for some objective function with the structure - min_s q_k(s) + ψ(x+s) s.t. ||s||_q⩽ Δ - for some Δ provided -Arguments ----------- -proxp : prox method for p-norm - takes in z (vector), a (λ||⋅||_p), p is norm for ψ I think -s0 : Vector{Float64,1} - Initial guess for the descent direction -projq : generic that projects onto ||⋅||_q⩽Δ norm ball -options : mutable structure p_params - - -Returns -------- -s : Vector{Float64,1} - Final value of Algorithm 6.1 descent direction -w : Vector{Float64,1} - relaxation variable of Algorithm 6.1 descent direction -""" -function prox_split_1w(proxp, s0, projq, options) - ε = options.optTol - max_iter = options.maxIter - ε_w = options.WoptTol - λ = options.λ - Δ = options.Δ - restart = options.restart - #ν_factor should have two values, as should ν - ν = options.ν - ξ = options.β^(-1) - η_factor = options.η_factor - ∇fk = options.∇fk - Bk = options.Bk - xk = options.xk - - if options.verbose == 0 - print_freq = Inf - elseif options.verbose == 1 - print_freq = round(max_iter / 10) - elseif options.verbose == 2 - print_freq = round(max_iter / 100) - else - print_freq = 1 - end - - s = zeros(size(s0)) - u = s0 + xk - w = zeros(size(s0)) - - err = 100 - w_err = norm(w - xk - s)^2 - k = 1 - - for i = 1:restart - while err > ε && k < max_iter && w_err > ε_w - #store previous values - s_ = s - u_ = u - w_ = w - # prox(z,α) = prox_lp(z, α, q) #x is what is changed, z is what is put in, p is norm, a is ξ - #I don't really like the Prox list that's currently in ProxLQ so I might make my own - - #u update with prox_(ξ*λ*||⋅||_p)(...) - u = proxp(u_ - ξ * (∇fk + Bk * (u_ - xk) - (w_ - u_ + xk) / η), ξ * λ) - - #update s - s = u - xk - - #w update - w = projq(s, Δ) - - w_err = norm(w - xk - s)^2 - err = norm(s_ - s) + norm(w_ - w) - k % print_freq == 0 && @printf( - "Iter: %4d, ||w-x-s||^2: %1.5e, s-err: %1.5e, w-err: %1.5e η: %1.5e, ξ: %1.5e \n", - k, - w_err, - norm(s_ - s), - norm(w_ - w), - η, - ξ - ) - - k = k + 1 - end - η *= η_factor - err = 100 - k = 0 - end - - return s, s_, w -end - -"""Solves descent direction s for some objective function with the structure - min_s q_k(s) + ψ(x+s) s.t. ||s||_q⩽ Δ - for some Δ provided -Arguments ----------- -proxp : prox method for p-norm - takes in z (vector), a (λ||⋅||_p), p is norm for ψ I think -s0 : Vector{Float64,1} - Initial guess for the descent direction -projq : generic that projects onto ||⋅||_q⩽Δ norm ball -options : mutable structure p_params - - -Returns -------- -s : Vector{Float64,1} - Final value of Algorithm 6.2 descent direction -w : Vector{Float64,1} - relaxation variable of Algorithm 6.2 descent direction -""" -function prox_split_2w(proxp, s0, projq, options) - ε = options.optTol - max_iter = options.maxIter - ε_w = options.WoptTol - λ = options.λ - Δ = options.Δ - restart = options.restart - #η_factor should have two values, as should η - η = options.η - ξ = η - η_factor = options.η_factor - ∇fk = options.∇fk - Bk = options.Bk - xk = options.xk - - print_freq = options.verbose - - u = s0 + xk - u_ = copy(u) - w1 = u #zeros(size(s0)) - w2 = u - xk #zeros(size(s0)) - - err = 100 - w1_err = norm(w1 - u)^2 - w2_err = norm(w2 - u + xk)^2 - s_feas = norm(s0, 1) - Δ - k = 1 - b_pt1 = Bk * xk - ∇fk - - for i = 1:restart - A = Bk + (2 / η) * I(size(Bk, 1)) - while err > ε && k < max_iter #&& (w1_err+w2_err)>ε_w - - #store previous values - u_ = u - w1_ = w1 - w2_ = w2 - # prox(z,α) = prox_lp(z, α, q) #x is what is changed, z is what is put in, p is norm, a is ξ - #I don't really like the Prox list that's currently in ProxLQ so I might make my own - - #u update with prox_(ξ*λ*||⋅||_p)(...) - # u = cg(A, b_pt1 + (w2+xk+w1)/η; maxiter=5) - u = fastcg(A, u_, b_pt1 + (w2 + xk + w1) / η; maxiter = 100) - - #update w1 - w1 = proxp(u, ξ * λ) - - #w2 update - w2 = projq(u - xk, Δ) - - w1_err = norm(w1 - u)^2 - w2_err = norm(w2 - u + xk)^2 - err = norm(u_ - u) + norm(w1_ - w1) + norm(w2_ - w2) - s_feas = norm(u - xk, 1) - Δ - k % print_freq == 0 && @printf( - "iter: %d, ||w1-u||²: %7.3e, ||w2-u+xk||²: %7.3e, err: %7.3e, η: %7.3e, s_feas: %7.3e, ||w1||_1: %7.3e, ||w2||_1: %7.3e, u-sum: %7.3e\n", - k, - w1_err, - w2_err, - err, - η, - s_feas, - norm(w1, 1), - norm(w2, 1), - ∇fk' * (u - xk) + 1 / 2 * (u - xk)' * Bk * (u - xk) - ) - - k = k + 1 - # η = η*η_factor - # ξ = η - end - η = η * η_factor - ξ = η - err = 100 - k = 1 - end - - return u - xk, u_ - xk, s_feas, err -end - -function fastcg(A, x, b; epsilon = 1.0f-5, maxiter = size(b, 1)) - r = b - A * x - p = r - rsold = r' * r - for i = 1:maxiter - Ap = A * p - alpha = rsold / (p' * Ap) - x = x + alpha * p - r = r - alpha * Ap - rsnew = r' * r - if sqrt(rsnew) < epsilon - break - end - p = r + (rsnew / rsold) * p - rsold = rsnew - end - return x -end