From 61ec6b11bf0d7a0ef3fc003f572cfc01acbf72b5 Mon Sep 17 00:00:00 2001 From: nantonel Date: Thu, 21 Feb 2019 18:07:45 +0100 Subject: [PATCH 1/3] rm blockarray stuff --- src/ProximalAlgorithms.jl | 2 +- .../AsymmetricForwardBackwardAdjoint.jl | 24 ++-- src/algorithms/DouglasRachford.jl | 18 +-- src/algorithms/ForwardBackward.jl | 94 ++++++-------- src/algorithms/PANOC.jl | 122 +++++++++--------- src/algorithms/ZeroFPR.jl | 99 +++++++------- src/template/Template.jl | 2 +- src/utilities/identity.jl | 4 +- src/utilities/tomove.jl | 10 ++ src/utilities/zero.jl | 8 +- test/REQUIRE | 3 +- test/runtests.jl | 7 +- test/test_lasso_small_split_f.jl | 16 +-- test/test_lasso_small_split_x.jl | 28 ++-- 14 files changed, 219 insertions(+), 218 deletions(-) create mode 100644 src/utilities/tomove.jl diff --git a/src/ProximalAlgorithms.jl b/src/ProximalAlgorithms.jl index 2be323c..0dfc51d 100644 --- a/src/ProximalAlgorithms.jl +++ b/src/ProximalAlgorithms.jl @@ -1,7 +1,6 @@ module ProximalAlgorithms using AbstractOperators -using AbstractOperators.BlockArrays using ProximalOperators using LinearAlgebra using Printf @@ -11,6 +10,7 @@ import Base: iterate include("utilities/identity.jl") include("utilities/zero.jl") include("utilities/conjugate.jl") +include("utilities/tomove.jl") abstract type ProximalAlgorithm{I,T} end diff --git a/src/algorithms/AsymmetricForwardBackwardAdjoint.jl b/src/algorithms/AsymmetricForwardBackwardAdjoint.jl index 471fd2c..49029bd 100644 --- a/src/algorithms/AsymmetricForwardBackwardAdjoint.jl +++ b/src/algorithms/AsymmetricForwardBackwardAdjoint.jl @@ -5,7 +5,7 @@ # [2] Condat. "A primal–dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms" Journal of Optimization Theory and Applications 158.2 (2013): 460-479. # [3] Vũ. "A splitting algorithm for dual monotone inclusions involving cocoercive operators"" Advances in Computational Mathematics, 38(3), pp.667-681. -struct AFBAIterator{I <: Integer, R <: Real, T1 <: BlockArray, T2 <: BlockArray} <: ProximalAlgorithm{I, Tuple{T1, T2}} +struct AFBAIterator{I <: Integer, R <: Real, T1 <: AbstractArray, T2 <: AbstractArray} <: ProximalAlgorithm{I, Tuple{T1, T2}} x::T1 y::T2 g @@ -37,19 +37,19 @@ end ################################################################################ # Constructor(s) -function AFBAIterator(x0::T1, y0::T2; g=IndFree(), h=IndFree(), f=IndFree(), l=IndZero(), L=Identity(blocksize(x0)), theta=1, mu=1, lam=1, betaQ=0, betaR=0, gamma1=-1, gamma2=-1, maxit::I=10000, tol::R=1e-5, verbose=1, verbose_freq = 100) where {I, R, T1, T2} - x = blockcopy(x0) - xbar = blockcopy(x0) - y = blockcopy(y0) - ybar = blockcopy(y0) - gradf = blockcopy(x0) - gradl = blockcopy(y0) - FPR_x = blockcopy(x0) +function AFBAIterator(x0::T1, y0::T2; g=IndFree(), h=IndFree(), f=IndFree(), l=IndZero(), L=Identity(size(x0)), theta=1, mu=1, lam=1, betaQ=0, betaR=0, gamma1=-1, gamma2=-1, maxit::I=10000, tol::R=1e-5, verbose=1, verbose_freq = 100) where {I, R, T1, T2} + x = copy(x0) + xbar = copy(x0) + y = copy(y0) + ybar = copy(y0) + gradf = copy(x0) + gradl = copy(y0) + FPR_x = copy(x0) FPR_x .= Inf - FPR_y = blockcopy(y0) + FPR_y = copy(y0) FPR_y .= Inf - temp_x = blockcopy(x0) - temp_y = blockcopy(y0) + temp_x = copy(x0) + temp_y = copy(y0) hconj = Conjugate(h) lconj = Conjugate(l) diff --git a/src/algorithms/DouglasRachford.jl b/src/algorithms/DouglasRachford.jl index b11b12d..6bf78a5 100644 --- a/src/algorithms/DouglasRachford.jl +++ b/src/algorithms/DouglasRachford.jl @@ -5,7 +5,7 @@ # [1] Eckstein, Bertsekas "On the Douglas-Rachford Splitting Method and the Proximal Point Algorithm for Maximal Monotone Operators*", Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989). # -struct DRSIterator{I <: Integer, R <: Real, T <: BlockArray{R}} <: ProximalAlgorithm{I,T} +struct DRSIterator{I <: Integer, R <: Real, T <: AbstractArray{R}} <: ProximalAlgorithm{I,T} x::T f g @@ -23,11 +23,11 @@ end ################################################################################ # Constructor(s) -function DRSIterator(x0::BlockArray{R}; f=Zero(), g=Zero(), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R} - y = blockcopy(x0) - r = blockcopy(x0) - z = blockcopy(x0) - FPR_x = blockcopy(x0) +function DRSIterator(x0::AbstractArray{R}; f=Zero(), g=Zero(), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R} + y = copy(x0) + r = copy(x0) + z = copy(x0) + FPR_x = copy(x0) FPR_x .= Inf return DRSIterator{I, R, typeof(x0)}(x0, f, g, gamma, maxit, tol, verbose, verbose_freq, y, r, z, FPR_x) end @@ -37,7 +37,7 @@ end maxit(sol::DRSIterator) = sol.maxit -converged(sol::DRSIterator, it) = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol +converged(sol::DRSIterator, it) = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol verbose(sol::DRSIterator) = sol.verbose > 0 verbose(sol::DRSIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0)) @@ -48,12 +48,12 @@ function display(sol::DRSIterator) end function display(sol::DRSIterator, it) - @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma) + @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma) end function Base.show(io::IO, sol::DRSIterator) println(io, "Douglas-Rachford Splitting" ) - println(io, "fpr : $(blockmaxabs(sol.FPR_x))") + println(io, "fpr : $(maximum(abs,sol.FPR_x))") print( io, "gamma : $(sol.gamma)") end diff --git a/src/algorithms/ForwardBackward.jl b/src/algorithms/ForwardBackward.jl index 6f5dd0d..e18eeed 100644 --- a/src/algorithms/ForwardBackward.jl +++ b/src/algorithms/ForwardBackward.jl @@ -50,31 +50,31 @@ end # Constructor function FBSIterator(x0::D; - fs::FS=Zero(), As::AS=Identity(blocksize(x0)), - fq::FQ=Zero(), Aq::AQ=Identity(blocksize(x0)), + fs::FS=Zero(), As::AS=Identity(size(x0)), + fq::FQ=Zero(), Aq::AQ=Identity(size(x0)), g::G=Zero(), gamma::R=-1.0, maxit::I=10000, tol::R=1e-4, adaptive::Bool=false, fast::Bool=false, verbose::I=1, verbose_freq::I = 100) where {I, R, D, FS, AS, FQ, AQ, G} - x = blockcopy(x0) - y = blockzeros(x0) - z = blockzeros(x0) - z_prev = blockzeros(x0) - FPR_x = blockzeros(x0) + x = copy(x0) + y = zero(x0) + z = zero(x0) + z_prev = zero(x0) + FPR_x = zero(x0) Aqx = Aq*x Asx = As*x - Aqz = blockzeros(Aqx) - Asz = blockzeros(Asx) - gradfq_Aqx = blockzeros(Aqx) - gradfs_Asx = blockzeros(Asx) - Aqt_gradfq_Aqx = blockzeros(x0) - Ast_gradfs_Asx = blockzeros(x0) - At_gradf_Ax = blockzeros(x0) - Aqz_prev = blockzeros(Aqx) - Asz_prev = blockzeros(Asx) - gradfq_Aqz = blockzeros(Aqx) - gradfs_Asz = blockzeros(Asx) - gradfq_Aqz_prev = blockzeros(Aqx) + Aqz = zero(Aqx) + Asz = zero(Asx) + gradfq_Aqx = zero(Aqx) + gradfs_Asx = zero(Asx) + Aqt_gradfq_Aqx = zero(x0) + Ast_gradfs_Asx = zero(x0) + At_gradf_Ax = zero(x0) + Aqz_prev = zero(Aqx) + Asz_prev = zero(Asx) + gradfq_Aqz = zero(Aqx) + gradfs_Asz = zero(Asx) + gradfq_Aqz_prev = zero(Aqx) CQ = typeof(Aqx) CS = typeof(Asx) FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}(x, @@ -96,7 +96,7 @@ end maxit(sol::FBSIterator{I}) where {I} = sol.maxit -converged(sol::FBSIterator{I, R, D}, it::I) where {I, R, D} = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol +converged(sol::FBSIterator{I, R, D}, it::I) where {I, R, D} = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol verbose(sol::FBSIterator) = sol.verbose > 0 verbose(sol::FBSIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0)) @@ -107,12 +107,12 @@ function display(sol::FBSIterator) end function display(sol::FBSIterator, it) - @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma) + @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma) end function Base.show(io::IO, sol::FBSIterator) println(io, (sol.fast ? "Fast " : "")*"Forward-Backward Splitting" ) - println(io, "fpr : $(blockmaxabs(sol.FPR_x))") + println(io, "fpr : $(maximum(abs,sol.FPR_x))") print( io, "gamma : $(sol.gamma)") end @@ -131,7 +131,7 @@ function initialize!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}) where sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx) mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx if sol.gamma <= 0.0 # estimate L in this case, and set gamma = 1/L @@ -146,15 +146,15 @@ function initialize!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}) where Asxeps = sol.As*xeps gradfs_Asxeps, = gradient(sol.fs, Asxeps) At_gradf_Axeps = sol.As'*gradfs_Asxeps .+ sol.Aq'*gradfq_Aqxeps - L = blockvecnorm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*blocklength(xeps))) + L = norm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*length(xeps))) sol.adaptive = true # in both cases set gamma = 1/L sol.gamma = 1.0/L end - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax prox!(sol.z, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z) + sol.FPR_x .= sol.x .- sol.z return sol.z end @@ -169,8 +169,8 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w if sol.adaptive for it_gam = 1:100 # TODO: replace/complement with lower bound on gamma - normFPR_x = blockvecnorm(sol.FPR_x) - uppbnd = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 + normFPR_x = norm(sol.FPR_x) + uppbnd = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 mul!(sol.Aqz, sol.Aq, sol.z) fq_Aqz = gradient!(sol.gradfq_Aqz, sol.fq, sol.Aqz) mul!(sol.Asz, sol.As, sol.z) @@ -178,9 +178,9 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w f_Az = fs_Asz + fq_Aqz if f_Az > uppbnd + 1e-6*abs(sol.f_Ax) sol.gamma = 0.5*sol.gamma - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x - sol.gamma .* sol.At_gradf_Ax prox!(sol.z, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z) + sol.FPR_x .= sol.x .- sol.z else break end @@ -198,7 +198,7 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w sol.gradfs_Asx, sol.gradfs_Asz = sol.gradfs_Asz, sol.gradfs_Asx mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx end else @@ -207,27 +207,19 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w extr = (sol.theta - 1.0)/theta1 sol.theta = theta1 # perform extrapolation - #sol.x .= sol.z .+ extr.*(sol.z .- sol.z_prev) - blockaxpy!(sol.x, sol.z, -1, sol.z_prev) - blockaxpy!(sol.x, sol.z, extr, sol.x) + sol.x .= sol.z .+ extr.*(sol.z .- sol.z_prev) sol.z, sol.z_prev = sol.z_prev, sol.z if sol.adaptive == true # extrapolate other extrapolable quantities - # diff_Aqx = (sol.Aqz .- sol.Aqz_prev) - # diff_gradfq_Aqx = (sol.gradfq_Aqz .- sol.gradfq_Aqz_prev) - blockaxpy!(sol.Aqx, sol.Aqz, -1.0, sol.Aqz_prev) - blockaxpy!(sol.gradfq_Aqx, sol.gradfq_Aqz, -1.0, sol.gradfq_Aqz_prev) - coeff_linear = real(blockvecdot(sol.gradfq_Aqz, sol.Aqx)) - coeff_quadr = real(blockvecdot(sol.Aqx, sol.gradfq_Aqx)) + sol.Aqx .= sol.Aqz .- sol.Aqz_prev + sol.gradfq_Aqx .= sol.gradfq_Aqz .- sol.gradfq_Aqz_prev + coeff_linear = real(dot(sol.gradfq_Aqz, sol.Aqx)) + coeff_quadr = real(dot(sol.Aqx, sol.gradfq_Aqx)) sol.fq_Aqx = fq_Aqz + extr*coeff_linear + 0.5*coeff_quadr*extr^2 - # sol.Aqx .= sol.Aqz .+ extr.*diff_Aqx - blockaxpy!(sol.Aqx, sol.Aqz, extr, sol.Aqx) - # sol.gradfq_Aqx .= sol.gradfq_Aqz .+ extr.*diff_gradfq_Aqx - blockaxpy!(sol.gradfq_Aqx, sol.gradfq_Aqz, extr, sol.gradfq_Aqx) - # sol.Asx .= sol.Asz .+ extr.*(sol.Asz .- sol.Asz_prev) - blockaxpy!(sol.Asx, sol.Asz, -1.0, sol.Asz_prev) - blockaxpy!(sol.Asx, sol.Asz, extr, sol.Asx ) + sol.Aqx .= sol.Aqz .+ extr.*(sol.Aqz .- sol.Aqz_prev) + sol.gradfq_Aqx .= sol.gradfq_Aqz .+ extr.*(sol.gradfq_Aqz .- sol.gradfq_Aqz_prev) + sol.Asx .= sol.Asz .+ extr .*(sol.Asz .- sol.Asz_prev) # store the z-quantities for future extrapolation sol.Aqz_prev, sol.Aqz = sol.Aqz, sol.Aqz_prev sol.gradfq_Aqz_prev, sol.gradfq_Aqz = sol.gradfq_Aqz, sol.gradfq_Aqz_prev @@ -237,7 +229,7 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) # TODO: we can probably save the MATVEC by Aq' in the next line mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx end end @@ -248,12 +240,12 @@ function iterate!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}, it::I) w sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx) mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx end - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax prox!(sol.z, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.z) + sol.FPR_x .= sol.x .- sol.z return sol.z diff --git a/src/algorithms/PANOC.jl b/src/algorithms/PANOC.jl index ddbf6ef..f2f1bbd 100644 --- a/src/algorithms/PANOC.jl +++ b/src/algorithms/PANOC.jl @@ -57,39 +57,39 @@ end # Constructor function PANOCIterator(x0::D; - fs::FS=Zero(), As::AS=Identity(blocksize(x0)), - fq::FQ=Zero(), Aq::AQ=Identity(blocksize(x0)), + fs::FS=Zero(), As::AS=Identity(size(x0)), + fq::FQ=Zero(), Aq::AQ=Identity(size(x0)), g::G=Zero(), gamma::R=-1.0, maxit::I=10000, tol::R=1e-4, adaptive::Bool=false, memory::I=10, verbose::I=1, verbose_freq::I=100, alpha::R=0.95, beta::R=0.5) where {I, R, D, FS, AS, FQ, AQ, G} - x = blockcopy(x0) - xbar = blockzeros(x0) - x_prev = blockzeros(x0) - xnew = blockzeros(x0) - xnewbar = blockzeros(x0) - y = blockzeros(x0) - FPR_x = blockzeros(x0) - FPR_x_prev = blockzeros(x0) - FPR_xnew = blockzeros(x0) + x = copy(x0) + xbar = zero(x0) + x_prev = zero(x0) + xnew = zero(x0) + xnewbar = zero(x0) + y = zero(x0) + FPR_x = zero(x0) + FPR_x_prev = zero(x0) + FPR_xnew = zero(x0) Aqx = Aq*x Asx = As*x - Aqxnew = blockzeros(Aqx) - Asxnew = blockzeros(Asx) - Aqd = blockzeros(Aqx) - Asd = blockzeros(Asx) - Aqfb = blockzeros(Aqx) - Asfb = blockzeros(Asx) - gradfq_Aqx = blockzeros(Aqx) - gradfs_Asx = blockzeros(Asx) - Aqxbar = blockzeros(Aqx) - Asxbar = blockzeros(Asx) - gradfq_Aqxbar = blockzeros(Aqx) - gradfs_Asxbar = blockzeros(Asx) - At_gradf_Ax = blockzeros(x0) - Aqt_gradfq_Aqx = blockzeros(x0) - Ast_gradfs_Asx = blockzeros(x0) - d = blockzeros(x0) + Aqxnew = zero(Aqx) + Asxnew = zero(Asx) + Aqd = zero(Aqx) + Asd = zero(Asx) + Aqfb = zero(Aqx) + Asfb = zero(Asx) + gradfq_Aqx = zero(Aqx) + gradfs_Asx = zero(Asx) + Aqxbar = zero(Aqx) + Asxbar = zero(Asx) + gradfq_Aqxbar = zero(Aqx) + gradfs_Asxbar = zero(Asx) + At_gradf_Ax = zero(x0) + Aqt_gradfq_Aqx = zero(x0) + Ast_gradfs_Asx = zero(x0) + d = zero(x0) H = LBFGS(x, memory) CQ = typeof(Aqx) CS = typeof(Asx) @@ -117,7 +117,7 @@ end maxit(sol::PANOCIterator{I}) where {I} = sol.maxit -converged(sol::PANOCIterator{I,R,D}, it::I) where {I,R,D}= it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol +converged(sol::PANOCIterator{I,R,D}, it::I) where {I,R,D}= it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol verbose(sol::PANOCIterator) = sol.verbose > 0 verbose(sol::PANOCIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0)) @@ -128,12 +128,12 @@ function display(sol::PANOCIterator) end function display(sol::PANOCIterator, it) - @printf("%6d | %7.4e | %7.4e | %7.4e | %7.4e | \n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma, sol.tau, sol.FBE_x) + @printf("%6d | %7.4e | %7.4e | %7.4e | %7.4e | \n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma, sol.tau, sol.FBE_x) end function Base.show(io::IO, sol::PANOCIterator) println(io, "PANOC" ) - println(io, "fpr : $(blockmaxabs(sol.FPR_x))") + println(io, "fpr : $(maximum(abs,sol.FPR_x))") println(io, "gamma : $(sol.gamma)") println(io, "tau : $(sol.tau)") print( io, "FBE : $(sol.FBE_x)") @@ -154,7 +154,7 @@ function initialize!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}) sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqx) mul!(sol.Asx, sol.As, sol.x) sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx) - blockaxpy!(sol.At_gradf_Ax, sol.As'*sol.gradfs_Asx, 1.0, sol.Aq'*sol.gradfq_Aqx) + sol.At_gradf_Ax .= sol.As'*sol.gradfs_Asx .+ sol.Aq'*sol.gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx if sol.gamma <= 0.0 # estimate L in this case, and set gamma = 1/L @@ -169,18 +169,18 @@ function initialize!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}) Asxeps = sol.As*xeps gradfs_Asxeps, = gradient(sol.fs, Asxeps) At_gradf_Axeps = sol.As'*gradfs_Asxeps .+ sol.Aq'*gradfq_Aqxeps - L = blockvecnorm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*blocklength(xeps))) + L = norm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*length(xeps))) sol.adaptive = true # in both cases set gamma = 1/L sol.gamma = sol.alpha/L end - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .-sol.gamma .* sol.At_gradf_Ax sol.g_xbar = prox!(sol.xbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.xbar) + sol.FPR_x .= sol.x .- sol.xbar - sol.normFPR_x = blockvecnorm(sol.FPR_x) - sol.FBE_x = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*sol.normFPR_x^2 + sol.g_xbar + sol.normFPR_x = norm(sol.FPR_x) + sol.FBE_x = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*sol.normFPR_x^2 + sol.g_xbar return sol.xbar @@ -199,14 +199,14 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it fs_Asxbar = gradient!(sol.gradfs_Asxbar, sol.fs, sol.Asxbar) f_Axbar = fs_Asxbar + fq_Aqxbar - uppbnd = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + + uppbnd = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*sol.normFPR_x^2 if f_Axbar > uppbnd + 1e-6*abs(sol.f_Ax) sol.gamma = 0.5*sol.gamma - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax sol.g_xbar = prox!(sol.xbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.xbar) - sol.normFPR_x = blockvecnorm(sol.FPR_x) + sol.FPR_x .= sol.x .- sol.xbar + sol.normFPR_x = norm(sol.FPR_x) else sol.FBE_x = uppbnd + sol.g_xbar break @@ -220,7 +220,7 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it mul!(sol.d, sol.H, ( x -> .-x ).(sol.FPR_x)) # TODO: not nice sol.FPR_x_prev, sol.FPR_x = sol.FPR_x, sol.FPR_x_prev - blockset!(sol.x_prev, sol.x) + sol.x_prev .= sol.x sigma = 0.5*sol.beta/sol.gamma*(1.0-sol.alpha) maxit_tau = 10 @@ -232,11 +232,11 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it mul!( sol.Asd, sol.As, sol.d) # xnew = x + tau*d - blockaxpy!(sol.xnew, sol.x, sol.tau, sol.d) + sol.xnew .= sol.x .+ sol.tau .* sol.d # Aq*xnew = Aq*x + tau*Aq*d - blockaxpy!(sol.Aqxnew, sol.Aqx, sol.tau, sol.Aqd) + sol.Aqxnew .= sol.Aqx .+ sol.tau .* sol.Aqd # As*xnew = As*x + tau*As*d - blockaxpy!(sol.Asxnew, sol.Asx, sol.tau, sol.Asd) + sol.Asxnew .= sol.Asx .+ sol.tau .* sol.Asd # calculate new FBE in xnew sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqxnew) @@ -245,18 +245,18 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) - blockaxpy!(sol.At_gradf_Ax, sol.Aqt_gradfq_Aqx, 1.0, sol.Ast_gradfs_Asx) + sol.At_gradf_Ax .= sol.Aqt_gradfq_Aqx .+ sol.Ast_gradfs_Asx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx # gradient step - blockaxpy!(sol.y, sol.xnew, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.xnew .- sol.gamma .* sol.At_gradf_Ax # prox step sol.g_xbar = prox!(sol.xnewbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_xnew, sol.xnew, -1.0, sol.xnewbar) - norm_FPRxnew = blockvecnorm(sol.FPR_xnew) + sol.FPR_xnew .= sol.xnew .- sol.xnewbar + norm_FPRxnew = norm(sol.FPR_xnew) - FBE_xnew = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_xnew)) + + FBE_xnew = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_xnew)) + 0.5/sol.gamma*norm_FPRxnew^2 + sol.g_xbar if FBE_xnew > sol.FBE_x - sigma*sol.normFPR_x^2 @@ -270,19 +270,19 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it sol.tau *= 0.5 # xnew = x + tau*d - blockaxpy!(sol.xnew, sol.x, sol.tau, sol.d) + sol.xnew .= sol.x .+ sol.tau .* sol.d # xnew = x + tau*d - (1-tau)*fb - blockaxpy!(sol.xnew, sol.xnew, sol.tau-1.0, sol.FPR_x_prev) + sol.xnew .= sol.xnew .+ (sol.tau-1.0) .* sol.FPR_x_prev # Aq*xnew = Aq*x + tau*Aq*d - blockaxpy!(sol.Aqxnew, sol.Aqx, sol.tau, sol.Aqd) + sol.Aqxnew .= sol.Aqx .+ sol.tau .* sol.Aqd # Aq*xnew = Aq*x + tau*Aq*d - (1-tau)*Aq*fb - blockaxpy!(sol.Aqxnew, sol.Aqxnew, sol.tau-1.0, sol.Aqfb) + sol.Aqxnew .= sol.Aqxnew .+ (sol.tau-1.0) .* sol.Aqfb # As*xnew = As*x + tau*As*d - blockaxpy!(sol.Asxnew, sol.Asx, sol.tau, sol.Asd) + sol.Asxnew .= sol.Asx .+ sol.tau .* sol.Asd # As*xnew = As*x + tau*As*d - (1-tau)*As*fb - blockaxpy!(sol.Asxnew, sol.Asxnew, sol.tau-1.0, sol.Asfb) + sol.Asxnew .= sol.Asxnew .+ (sol.tau-1.0) .* sol.Asfb # calculate new FBE in xnew sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqxnew) @@ -291,18 +291,18 @@ function iterate!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, it mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) - blockaxpy!(sol.At_gradf_Ax, sol.Aqt_gradfq_Aqx, 1.0, sol.Ast_gradfs_Asx) + sol.At_gradf_Ax .= sol.Aqt_gradfq_Aqx .+ sol.Ast_gradfs_Asx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx # gradient step - blockaxpy!(sol.y, sol.xnew, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.xnew .- sol.gamma .* sol.At_gradf_Ax # prox step sol.g_xbar = prox!(sol.xnewbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_xnew, sol.xnew, -1.0, sol.xnewbar) - norm_FPRxnew = blockvecnorm(sol.FPR_xnew) + sol.FPR_xnew .= sol.xnew .- sol.xnewbar + norm_FPRxnew = norm(sol.FPR_xnew) - FBE_xnew = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_xnew)) + + FBE_xnew = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_xnew)) + 0.5/sol.gamma*norm_FPRxnew^2 + sol.g_xbar if FBE_xnew <= sol.FBE_x - sigma*sol.normFPR_x^2 @@ -330,6 +330,6 @@ end function PANOC(x0; kwargs...) sol = PANOCIterator(x0; kwargs...) it, point = run!(sol) - blockset!(x0, point) + x0 .= point return (it, point, sol) end diff --git a/src/algorithms/ZeroFPR.jl b/src/algorithms/ZeroFPR.jl index 3512078..4dc0e17 100644 --- a/src/algorithms/ZeroFPR.jl +++ b/src/algorithms/ZeroFPR.jl @@ -51,35 +51,35 @@ end # Constructor function ZeroFPRIterator(x0::D; - fs::FS=Zero(), As::AS=Identity(blocksize(x0)), - fq::FQ=Zero(), Aq::AQ=Identity(blocksize(x0)), + fs::FS=Zero(), As::AS=Identity(size(x0)), + fq::FQ=Zero(), Aq::AQ=Identity(size(x0)), g::G=Zero(), gamma::R=-1.0, maxit::I=10000, tol::R=1e-4, adaptive::Bool=false, memory::I=10, verbose::I=1, verbose_freq::I=100, alpha::R=0.95, beta::R=0.5) where {I, R, D, FS, AS, FQ, AQ, G} - x = blockcopy(x0) - y = blockzeros(x0) - xbar = blockzeros(x0) - xbarbar = blockzeros(x0) - xnew = blockzeros(x0) - xnewbar = blockzeros(x0) - xbar = blockzeros(x0) - FPR_x = blockzeros(x0) - FPR_xbar_prev = blockzeros(x0) - FPR_xbar = blockzeros(x0) + x = copy(x0) + y = zero(x0) + xbar = zero(x0) + xbarbar = zero(x0) + xnew = zero(x0) + xnewbar = zero(x0) + xbar = zero(x0) + FPR_x = zero(x0) + FPR_xbar_prev = zero(x0) + FPR_xbar = zero(x0) Aqx = Aq*x Asx = As*x - Aqxbar = blockzeros(Aqx) - Asxbar = blockzeros(Asx) - Aqxnew = blockzeros(Aqx) - Asxnew = blockzeros(Asx) - Aqd = blockzeros(Aqx) - Asd = blockzeros(Asx) - gradfq_Aqx = blockzeros(Aqx) - gradfs_Asx = blockzeros(Asx) - At_gradf_Ax = blockzeros(x0) - Ast_gradfs_Asx = blockzeros(x0) - Aqt_gradfq_Aqx = blockzeros(x0) - d = blockzeros(x0) + Aqxbar = zero(Aqx) + Asxbar = zero(Asx) + Aqxnew = zero(Aqx) + Asxnew = zero(Asx) + Aqd = zero(Aqx) + Asd = zero(Asx) + gradfq_Aqx = zero(Aqx) + gradfs_Asx = zero(Asx) + At_gradf_Ax = zero(x0) + Ast_gradfs_Asx = zero(x0) + Aqt_gradfq_Aqx = zero(x0) + d = zero(x0) H = LBFGS(x, memory) CQ = typeof(Aqx) CS = typeof(Asx) @@ -107,7 +107,7 @@ end maxit(sol::ZeroFPRIterator{I}) where {I} = sol.maxit -converged(sol::ZeroFPRIterator{I,R,D}, it::I) where {I,R,D} = it > 0 && blockmaxabs(sol.FPR_x)/sol.gamma <= sol.tol +converged(sol::ZeroFPRIterator{I,R,D}, it::I) where {I,R,D} = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol verbose(sol::ZeroFPRIterator) = sol.verbose > 0 verbose(sol::ZeroFPRIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0)) @@ -117,12 +117,12 @@ function display(sol::ZeroFPRIterator) @printf("------|------------|------------|------------|------------|\n") end function display(sol::ZeroFPRIterator, it) - @printf("%6d | %7.4e | %7.4e | %7.4e | %7.4e | \n", it, sol.gamma, blockmaxabs(sol.FPR_x)/sol.gamma, sol.tau, sol.FBE_x) + @printf("%6d | %7.4e | %7.4e | %7.4e | %7.4e | \n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma, sol.tau, sol.FBE_x) end function Base.show(io::IO, sol::ZeroFPRIterator) println(io, "ZeroFPR" ) - println(io, "fpr : $(blockmaxabs(sol.FPR_x))") + println(io, "fpr : $(maximum(abs,sol.FPR_x))") println(io, "gamma : $(sol.gamma)") println(io, "tau : $(sol.tau)") print( io, "FBE : $(sol.FBE_x)") @@ -143,7 +143,7 @@ function initialize!(sol::ZeroFPRIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH sol.fq_Aqx = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqx) mul!(sol.Asx, sol.As, sol.x) sol.fs_Asx = gradient!(sol.gradfs_Asx, sol.fs, sol.Asx) - blockaxpy!(sol.At_gradf_Ax, sol.As'*sol.gradfs_Asx, 1.0, sol.Aq'*sol.gradfq_Aqx) + sol.At_gradf_Ax .= sol.As'*sol.gradfs_Asx .+ sol.Aq'*sol.gradfq_Aqx sol.f_Ax = sol.fs_Asx + sol.fq_Aqx if sol.gamma <= 0.0 # estimate L in this case, and set gamma = 1/L @@ -151,22 +151,21 @@ function initialize!(sol::ZeroFPRIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH # 1) if adaptive = false and only fq is present then L is "accurate" # 2) otherwise L is "inaccurate" and set adaptive = true # TODO: implement case 1), now 2) is always performed - # xeps = sol.x .+ sqrt(eps()) - xeps = (x -> x .+ sqrt(eps())).(sol.x) + xeps = sol.x .+ sqrt(eps()) Aqxeps = sol.Aq*xeps gradfq_Aqxeps, = gradient(sol.fq, Aqxeps) Asxeps = sol.As*xeps gradfs_Asxeps, = gradient(sol.fs, Asxeps) At_gradf_Axeps = sol.As'*gradfs_Asxeps .+ sol.Aq'*gradfq_Aqxeps - L = blockvecnorm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*blocklength(xeps))) + L = norm(sol.At_gradf_Ax .- At_gradf_Axeps)/(sqrt(eps()*length(xeps))) sol.adaptive = true # in both cases set gamma = 1/L sol.gamma = sol.alpha/L end - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax sol.g_xbar = prox!(sol.xbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.xbar) + sol.FPR_x .= sol.x .- sol.xbar return sol.xbar @@ -186,13 +185,13 @@ function iterate!(sol::ZeroFPRIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, if sol.adaptive for it_gam = 1:100 # TODO: replace/complement with lower bound on gamma - normFPR_x = blockvecnorm(sol.FPR_x) - uppbnd = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 + normFPR_x = norm(sol.FPR_x) + uppbnd = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 if f_Axbar > uppbnd + 1e-6*abs(sol.f_Ax) sol.gamma = 0.5*sol.gamma - blockaxpy!(sol.y, sol.x, -sol.gamma, sol.At_gradf_Ax) + sol.y .= sol.x .- sol.gamma .* sol.At_gradf_Ax sol.xbar, sol.g_xbar = prox(sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.x, -1.0, sol.xbar) + sol.FPR_x .= sol.x .- sol.xbar else break end @@ -206,16 +205,16 @@ function iterate!(sol::ZeroFPRIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, # Compute value of FBE at x - normFPR_x = blockvecnorm(sol.FPR_x) - FBE_x = sol.f_Ax - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 + sol.g_xbar + normFPR_x = norm(sol.FPR_x) + FBE_x = sol.f_Ax - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_x^2 + sol.g_xbar # Compute search direction mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) - blockaxpy!(sol.y, sol.xbar, -sol.gamma, sol.At_gradf_Ax) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx + sol.y .= sol.xbar .- sol.gamma .* sol.At_gradf_Ax g_xbarbar = prox!(sol.xbarbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_xbar, sol.xbar, -1.0, sol.xbarbar) + sol.FPR_xbar .= sol.xbar .- sol.xbarbar if it > 1 update!(sol.H, sol.xbar, sol.xnewbar, sol.FPR_xbar, sol.FPR_xbar_prev) @@ -237,21 +236,21 @@ function iterate!(sol::ZeroFPRIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}, maxit_tau = 10 for it_tau = 1:maxit_tau # TODO: replace/complement with lower bound on tau - blockaxpy!(sol.xnew, sol.xbar, sol.tau, sol.d) - blockaxpy!(sol.Asxnew, sol.Asxbar, sol.tau, sol.Asd) - blockaxpy!(sol.Aqxnew, sol.Aqxbar, sol.tau, sol.Aqd) + sol.xnew .= sol.xbar .+ sol.tau .* sol.d + sol.Asxnew .= sol.Asxbar .+ sol.tau .* sol.Asd + sol.Aqxnew .= sol.Aqxbar .+ sol.tau .* sol.Aqd fs_Asxnew = gradient!(sol.gradfs_Asx, sol.fs, sol.Asxnew) # TODO: can precompute most of next line before the iteration fq_Aqxnew = gradient!(sol.gradfq_Aqx, sol.fq, sol.Aqxnew) f_Axnew = fs_Asxnew + fq_Aqxnew mul!(sol.Ast_gradfs_Asx, sol.As', sol.gradfs_Asx) mul!(sol.Aqt_gradfq_Aqx, sol.Aq', sol.gradfq_Aqx) - blockaxpy!(sol.At_gradf_Ax, sol.Ast_gradfs_Asx, 1.0, sol.Aqt_gradfq_Aqx) - blockaxpy!(sol.y, sol.xnew, -sol.gamma, sol.At_gradf_Ax) + sol.At_gradf_Ax .= sol.Ast_gradfs_Asx .+ sol.Aqt_gradfq_Aqx + sol.y .= sol.xnew .- sol.gamma .* sol.At_gradf_Ax g_xnewbar = prox!(sol.xnewbar, sol.g, sol.y, sol.gamma) - blockaxpy!(sol.FPR_x, sol.xnew, -1.0, sol.xnewbar) - normFPR_xnew = blockvecnorm(sol.FPR_x) - FBE_xnew = f_Axnew - real(blockvecdot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_xnew^2 + g_xnewbar + sol.FPR_x .= sol.xnew .- sol.xnewbar + normFPR_xnew = norm(sol.FPR_x) + FBE_xnew = f_Axnew - real(dot(sol.At_gradf_Ax, sol.FPR_x)) + 0.5/sol.gamma*normFPR_xnew^2 + g_xnewbar if FBE_xnew <= FBE_x - sigma*normFPR_x^2 break end diff --git a/src/template/Template.jl b/src/template/Template.jl index c4a3c2e..a5f1839 100644 --- a/src/template/Template.jl +++ b/src/template/Template.jl @@ -1,7 +1,7 @@ ################################################################################ # Template iterator -struct TemplateIterator{I <: Integer, T <: BlockArray} <: ProximalAlgorithm{I,T} +struct TemplateIterator{I <: Integer, T <: AbstractArray} <: ProximalAlgorithm{I,T} x::T maxit::I # Put here problem description, other parameters, method memory, diff --git a/src/utilities/identity.jl b/src/utilities/identity.jl index e6e16af..dc1c91c 100644 --- a/src/utilities/identity.jl +++ b/src/utilities/identity.jl @@ -8,11 +8,11 @@ end size(A::Identity) = (A.size, A.size) size(A::Identity, dim::Integer) = A.size -mul!(y, A::Identity, x) = blockcopy!(y, x) +mul!(y, A::Identity, x) = y .= x adjoint(A::Identity) = A -(*)(A::Identity, x) = blockcopy(x) +(*)(A::Identity, x) = copy(x) opnorm(A::Identity) = 1.0 opnorm(A::Identity, p::Real) = 1.0 diff --git a/src/utilities/tomove.jl b/src/utilities/tomove.jl new file mode 100644 index 0000000..65380d3 --- /dev/null +++ b/src/utilities/tomove.jl @@ -0,0 +1,10 @@ +####### TO MOVE +using RecursiveArrayTools +function prox(f::ProximableFunction, x::ArrayPartition, args...) + y, fy = prox(f, x.x, args...) + return ArrayPartition(y), fy +end + +prox!(y::ArrayPartition, f::ProximableFunction, x::ArrayPartition, args...) = +prox!(y.x, f, x.x, args...) +############### diff --git a/src/utilities/zero.jl b/src/utilities/zero.jl index 5b52e42..748b75f 100644 --- a/src/utilities/zero.jl +++ b/src/utilities/zero.jl @@ -5,21 +5,21 @@ struct Zero end (f::Zero)(x) = 0.0 function gradient!(y, f::Zero, x) - blockset!(y, 0.0) + fill!(y,0.0) return 0.0 end function gradient(f::Zero, x) - y = blockzeros(blocksize(x)) + y = zeros(size(x)) return y, 0.0 end function prox!(y, f::Zero, x, gamma) - blockcopy!(y, x) + y .= x return 0.0 end function prox(f::Zero, x, gamma) - y = blockcopy(x) + y = copy(x) return y, 0.0 end diff --git a/test/REQUIRE b/test/REQUIRE index d54f1de..1aeb5ea 100644 --- a/test/REQUIRE +++ b/test/REQUIRE @@ -1,3 +1,4 @@ -julia 0.6 +julia 1.0 AbstractOperators 0.0.3 ProximalOperators 0.6.0 +RecursiveArrayTools 0.18 diff --git a/test/runtests.jl b/test/runtests.jl index b64fdb9..2ab2fea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,6 @@ using ProximalOperators using ProximalAlgorithms using AbstractOperators -using AbstractOperators.BlockArrays using LinearAlgebra using SparseArrays using FFTW @@ -9,10 +8,10 @@ using Random using Test using DelimitedFiles + @testset "ProximalAlgorithms" begin @testset "Utilities" begin -# include("test_block.jl") # remove? include("test_conjugate.jl") end @@ -20,8 +19,8 @@ end include("test_template.jl") include("test_lasso_small.jl") include("test_lasso_harder.jl") - include("test_lasso_small_split_x.jl") - include("test_lasso_small_split_f.jl") +# include("test_lasso_small_split_x.jl") +# include("test_lasso_small_split_f.jl") include("test_l1logreg_small.jl") include("test_afba.jl") include("test_afba_LP.jl") diff --git a/test/test_lasso_small_split_f.jl b/test/test_lasso_small_split_f.jl index e2da711..fb80059 100644 --- a/test/test_lasso_small_split_f.jl +++ b/test/test_lasso_small_split_f.jl @@ -20,7 +20,7 @@ x_star = [-3.877278911564627e-01, 0, 0, 2.174149659863943e-02, 6.168435374149660 # Nonfast/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2) @test norm(x .- x_star, Inf) <= 1e-4 @test it < 150 @@ -28,7 +28,7 @@ println(sol) # Nonfast/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true) @test norm(x .- x_star, Inf) <= 1e-4 @test it < 300 @@ -36,7 +36,7 @@ println(sol) # Fast/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2, fast=true) @test norm(x .- x_star, Inf) <= 1e-4 @test it < 100 @@ -44,7 +44,7 @@ println(sol) # Fast/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true, fast=true) @test norm(x .- x_star, Inf) <= 1e-4 @test it < 200 @@ -52,7 +52,7 @@ println(sol) # ZeroFPR/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2) @test norm(x - x_star, Inf) <= 1e-4 @test it < 15 @@ -60,14 +60,14 @@ println(sol) # ZeroFPR/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, adaptive=true) @test norm(x - x_star, Inf) <= 1e-4 @test it < 20 # PANOC/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2) @test norm(x - x_star, Inf) <= 1e-4 @test it < 20 @@ -75,7 +75,7 @@ println(sol) # PANOC/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, adaptive=true) @test norm(x - x_star, Inf) <= 1e-4 @test it < 20 diff --git a/test/test_lasso_small_split_x.jl b/test/test_lasso_small_split_x.jl index 3759609..2405b0a 100644 --- a/test/test_lasso_small_split_x.jl +++ b/test/test_lasso_small_split_x.jl @@ -1,4 +1,4 @@ -using AbstractOperators +using AbstractOperators, RecursiveArrayTools using ProximalOperators A1 = [ 1.0 -2.0 3.0; @@ -23,52 +23,52 @@ f = Translate(SqrNormL2(), -b) lam = 0.1*norm(A'*b, Inf) g = SeparableSum(NormL1(lam), NormL1(lam)) -x_star = ([-3.877278911564627e-01, 0, 0], [2.174149659863943e-02, 6.168435374149660e-01]) +x_star = ArrayPartition([-3.877278911564627e-01, 0, 0], [2.174149659863943e-02, 6.168435374149660e-01]) # Nonfast/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs,x .- x_star) <= 1e-4 @test it < 150 println(sol) # Nonfast/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs,x .- x_star) <= 1e-4 @test it < 300 println(sol) # Fast/Nonadaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, gamma=1.0/opnorm(A)^2, fast=true) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs, x .- x_star) <= 1e-4 @test it < 100 println(sol) # Fast/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.FBS(x0; fq=f, Aq=opA, g=g, adaptive=true, fast=true) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs,x .- x_star) <= 1e-4 @test it < 200 println(sol) # ZeroFPR/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.ZeroFPR(x0; fq=f, Aq=opA, g=g, adaptive=true) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs,x .- x_star) <= 1e-4 @test it < 15 println(sol) # PANOC/Adaptive -x0 = ProximalAlgorithms.blockzeros(x_star) +x0 = zero(x_star) @time it, x, sol = ProximalAlgorithms.PANOC(x0; fq=f, Aq=opA, g=g, adaptive=true) -@test ProximalAlgorithms.blockmaxabs(x .- x_star) <= 1e-4 +@test maximum(abs,x .- x_star) <= 1e-4 @test it < 20 println(sol) From 1c710269f2aa523fa0a85b5cc47555f5c45a1567 Mon Sep 17 00:00:00 2001 From: nantonel Date: Fri, 22 Feb 2019 09:48:44 +0100 Subject: [PATCH 2/3] all test passing --- src/algorithms/ForwardBackward.jl | 3 +-- src/algorithms/PANOC.jl | 3 +-- src/utilities/tomove.jl | 10 ++++++++++ src/utilities/zero.jl | 2 +- test/runtests.jl | 4 ++-- 5 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/algorithms/ForwardBackward.jl b/src/algorithms/ForwardBackward.jl index e18eeed..443ce13 100644 --- a/src/algorithms/ForwardBackward.jl +++ b/src/algorithms/ForwardBackward.jl @@ -139,8 +139,7 @@ function initialize!(sol::FBSIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G}) where # 1) if adaptive = false and only fq is present then L is "accurate" # 2) otherwise L is "inaccurate" and set adaptive = true # TODO: implement case 1), now 2) is always performed - # xeps = sol.x .+ sqrt(eps()) - xeps = (x -> x .+ sqrt(eps())).(sol.x) + xeps = sol.x .+ sqrt(eps()) Aqxeps = sol.Aq*xeps gradfq_Aqxeps, = gradient(sol.fq, Aqxeps) Asxeps = sol.As*xeps diff --git a/src/algorithms/PANOC.jl b/src/algorithms/PANOC.jl index f2f1bbd..b08e44c 100644 --- a/src/algorithms/PANOC.jl +++ b/src/algorithms/PANOC.jl @@ -162,8 +162,7 @@ function initialize!(sol::PANOCIterator{I, R, D, CS, FS, AS, CQ, FQ, AQ, G, HH}) # 1) if adaptive = false and only fq is present then L is "accurate" # 2) otherwise L is "inaccurate" and set adaptive = true # TODO: implement case 1), now 2) is always performed - # xeps = sol.x .+ sqrt(eps()) - xeps = (x -> x .+ sqrt(eps())).(sol.x) + xeps = sol.x .+ sqrt(eps()) Aqxeps = sol.Aq*xeps gradfq_Aqxeps, = gradient(sol.fq, Aqxeps) Asxeps = sol.As*xeps diff --git a/src/utilities/tomove.jl b/src/utilities/tomove.jl index 65380d3..e33aafa 100644 --- a/src/utilities/tomove.jl +++ b/src/utilities/tomove.jl @@ -1,4 +1,6 @@ ####### TO MOVE +import ProximalOperators: prox, prox!, gradient, gradient! + using RecursiveArrayTools function prox(f::ProximableFunction, x::ArrayPartition, args...) y, fy = prox(f, x.x, args...) @@ -7,4 +9,12 @@ end prox!(y::ArrayPartition, f::ProximableFunction, x::ArrayPartition, args...) = prox!(y.x, f, x.x, args...) + +function gradient(f::ProximableFunction, x::ArrayPartition) + gradx, fy = gradient(f, x.x) + return ArrayPartition(gradx), fy +end + +gradient!(y::ArrayPartition, f::ProximableFunction, x::ArrayPartition) = +gradient!(y.x, f, x.x) ############### diff --git a/src/utilities/zero.jl b/src/utilities/zero.jl index 748b75f..ccf25e1 100644 --- a/src/utilities/zero.jl +++ b/src/utilities/zero.jl @@ -10,7 +10,7 @@ function gradient!(y, f::Zero, x) end function gradient(f::Zero, x) - y = zeros(size(x)) + y = zero(x) return y, 0.0 end diff --git a/test/runtests.jl b/test/runtests.jl index 2ab2fea..f231de7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,8 +19,8 @@ end include("test_template.jl") include("test_lasso_small.jl") include("test_lasso_harder.jl") -# include("test_lasso_small_split_x.jl") -# include("test_lasso_small_split_f.jl") + include("test_lasso_small_split_x.jl") + include("test_lasso_small_split_f.jl") include("test_l1logreg_small.jl") include("test_afba.jl") include("test_afba_LP.jl") From 29629f31ccd179ad581bdf514f6969e5e28d1b92 Mon Sep 17 00:00:00 2001 From: nantonel Date: Fri, 22 Feb 2019 10:05:20 +0100 Subject: [PATCH 3/3] dropped 0.7 Travis and appveyor --- .travis.yml | 2 +- appveyor.yml | 3 +-- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.travis.yml b/.travis.yml index 0cc4353..8fb79df 100644 --- a/.travis.yml +++ b/.travis.yml @@ -4,8 +4,8 @@ os: - linux - osx julia: - - 0.7 - 1.0 + - 1.1 - nightly matrix: allow_failures: diff --git a/appveyor.yml b/appveyor.yml index 77bae98..27b4f11 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -1,7 +1,7 @@ environment: matrix: - - julia_version: 0.7 - julia_version: 1 + - julia_version: 1.1 - julia_version: nightly platform: @@ -10,7 +10,6 @@ platform: matrix: allow_failures: - - julia_version: 1 - julia_version: nightly branches: