update solvers to follow ProximalAlgorithms

docs/src/solvers.md

@@ -9,7 +9,7 @@
 !!! note "Problem warm-starting"
 
     By default *warm-starting* is always enabled.
-    For example, if two problems that utilize the same variables are solved consecutively,
+    For example, if two problems that involve the same variables are solved consecutively,
     the second one will be automatically warm-started by the solution of the first one.
     That is because the variables are always linked to their respective data vectors.
     If one wants to avoid this, the optimization variables needs to be manually re-initialized
@@ -18,20 +18,19 @@
 
 ## Specifying solver and options
 
-As shown above it is possible to choose the type of algorithm and specify its options by creating a `Solver` object.
-Currently, the following algorithms are supported:
+You can pick the algorithm to use as `Solver` object from the
+[`ProximalAlgorithms.jl`](https://github.com/kul-forbes/ProximalAlgorithms.jl))
+package. Currently, the following algorithms are supported:
 
-* *Proximal Gradient (PG)* [[1]](http://www.mit.edu/~dimitrib/PTseng/papers/apgm.pdf), [[2]](http://epubs.siam.org/doi/abs/10.1137/080716542)
-* *Fast Proximal Gradient (FPG)* [[1]](http://www.mit.edu/~dimitrib/PTseng/papers/apgm.pdf), [[2]](http://epubs.siam.org/doi/abs/10.1137/080716542)
-* *ZeroFPR* [[3]](https://arxiv.org/abs/1606.06256)
-* *PANOC* [[4]](https://doi.org/10.1109/CDC.2017.8263933)
-
-```@docs
-PG
-FPG
-ZeroFPR
-PANOC
-```
+* `ProximalAlgorithms.ForwardBackward`, also known as *proximal gradient*
+method [[1]](http://www.mit.edu/~dimitrib/PTseng/papers/apgm.pdf), [[2]](http://epubs.siam.org/doi/abs/10.1137/080716542). Nesterov acceleration can be enabled, which significantly
+improves its performance for convex problems.
+* `ProximalAlgorithms.ZeroFPR`, a Newton-type forward-backward algorithm,
+proposed in [[3]](https://arxiv.org/abs/1606.06256), using L-BFGS
+directions to accelerate convergence.
+* `ProximalAlgorithms.PANOC`, another Newton-type forward-backward algorithm,
+proposed in [[4]](https://doi.org/10.1109/CDC.2017.8263933), also using
+L-BFGS directions.
 
 ## Build and solve
 
@@ -43,18 +42,8 @@ problem
 solve
 ```
 
-It is important to stress out that the `Solver` objects created using
-the functions above ([`PG`](@ref), [`FPG`](@ref), etc.)
-specify only the type of algorithm to be used together with its options.
-The actual solver
-(namely the one of [`ProximalAlgorithms.jl`](https://github.com/kul-forbes/ProximalAlgorithms.jl))
-is constructed altogether with the problem formulation.
-The problem parsing procedure can be separated from the solver application using the functions [`build`](@ref) and [`solve!`](@ref).
-
-```@docs
-build
-solve!
-```
+Once again, the `Solver` objects is to be picked from
+[`ProximalAlgorithms.jl`](https://github.com/kul-forbes/ProximalAlgorithms.jl)).
 
 ## References
 

src/StructuredOptimization.jl

@@ -9,7 +9,18 @@ using ProximalOperators
 using ProximalAlgorithms
 
 include("syntax/syntax.jl")
-include("calculus/precomposeNonlinear.jl") #TODO move to ProximalOperators?
-include("solvers/solvers.jl")
+include("calculus/precomposeNonlinear.jl") # TODO move to ProximalOperators?
+include("arraypartition.jl") # TODO move to ProximalOperators?
+
+# problem parsing
+include("solvers/terms_extract.jl")
+include("solvers/terms_properties.jl")
+include("solvers/terms_splitting.jl")
+
+# solver calls
+include("solvers/solvers_options.jl")
+include("solvers/build_solve.jl")
+include("solvers/minimize.jl")
+
 
 end

src/arraypartition.jl

@@ -0,0 +1,36 @@
+import ProximalOperators
+import RecursiveArrayTools
+
+@inline function ProximalOperators.prox(
+    h::ProximalOperators.ProximableFunction,
+    x::RecursiveArrayTools.ArrayPartition,
+    gamma...
+)
+    # unwrap
+    y, fy = ProximalOperators.prox(h, x.x, gamma...)
+    # wrap
+    return RecursiveArrayTools.ArrayPartition(y), fy
+end
+
+@inline function ProximalOperators.gradient(
+    h::ProximalOperators.ProximableFunction,
+    x::RecursiveArrayTools.ArrayPartition
+)
+    # unwrap
+    grad, fx = ProximalOperators.gradient(h, x.x)
+    # wrap
+    return RecursiveArrayTools.ArrayPartition(grad), fx
+end
+
+@inline ProximalOperators.prox!(
+    y::RecursiveArrayTools.ArrayPartition,
+    h::ProximalOperators.ProximableFunction,
+    x::RecursiveArrayTools.ArrayPartition,
+    gamma...
+) = ProximalOperators.prox!(y.x, h, x.x, gamma...)
+
+@inline ProximalOperators.gradient!(
+    y::RecursiveArrayTools.ArrayPartition,
+    h::ProximalOperators.ProximableFunction,
+    x::RecursiveArrayTools.ArrayPartition
+) = ProximalOperators.gradient!(y.x, h, x.x)

src/solvers/build_solve.jl

@@ -1,11 +1,12 @@
 export build
 
 """
-`build(terms::Tuple, solver_opt::ForwardBackwardSolver)`
+    parse_problem(terms::Tuple, solver::ForwardBackwardSolver)
 
-Takes as input a tuple containing the terms defining the problem and the solver options.
+Takes as input a tuple containing the terms defining the problem and the solver.
 
-Returns a tuple containing the optimization variables and the built solver.
+Returns a tuple containing the optimization variables and the problem terms
+to be fed into the solver.
 
 # Example
 
@@ -18,82 +19,37 @@ julia> A, b = randn(10,4), randn(10);
 julia> p = problem( ls(A*x - b ) , norm(x) <= 1 );
 
 julia> build(p, PG());
-
 ```
-
 """
-function build(terms::Tuple, solver::ForwardBackwardSolver)
+function parse_problem(terms::Tuple, solver::T) where T <: ForwardBackwardSolver
   x = extract_variables(terms)
   # Separate smooth and nonsmooth
   smooth, nonsmooth = split_smooth(terms)
-  # Separate quadratic and nonquadratic
-  quadratic, smooth = split_quadratic(smooth)
-  kwargs = Array{Any, 1}()
   if is_proximable(nonsmooth)
     g = extract_proximable(x, nonsmooth)
-    append!(kwargs, [(:g, g)])
-    if !isempty(quadratic)
-      fq = extract_functions(quadratic)
-      Aq = extract_operators(x, quadratic)
-      append!(kwargs, [(:fq, fq)])
-      append!(kwargs, [(:Aq, Aq)])
-    end
+    kwargs = Dict{Symbol, Any}(:g => g)
     if !isempty(smooth)
       if is_linear(smooth)
-        fs = extract_functions(smooth)
-        As = extract_operators(x, smooth)
-        append!(kwargs, [(:As, As)])
-      else
-        fs = extract_functions_nodisp(smooth)
-        As = extract_affines(x, smooth)
-        fs = PrecomposeNonlinear(fs, As)
+        f = extract_functions(smooth)
+        A = extract_operators(x, smooth)
+        kwargs[:A] = A
+      else  # ??
+        f = extract_functions_nodisp(smooth)
+        A = extract_affines(x, smooth)
+        f = PrecomposeNonlinear(f, A)
       end
-      append!(kwargs, [(:fs, fs)])
+      kwargs[:f] = f
     end
-    return build_iterator(x, solver; kwargs...)
+    return (x, kwargs)
   end
-  error("Sorry, I cannot solve this problem")
-end
-
-################################################################################
-export solve!
-
-"""
-`solve!( x_solver )`
-
-Takes as input a tuple containing the optimization variables and the built solver.
-
-Solves the problem returning a tuple containing the iterations taken and the build solver.
-
-# Example
-
-```julia
-julia> x = Variable(4)
-Variable(Float64, (4,))
-
-julia> A, b = randn(10,4), randn(10);
-
-julia> p = problem( ls(A*x - b ) , norm(x) <= 1 );
-
-julia> x_solver = build(p, PG(verbose = 0));
-
-julia> solve!(x_solver);
-
-```
-
-"""
-function solve!(x_and_iter::Tuple{Tuple{Vararg{Variable}}, ProximalAlgorithms.ProximalAlgorithm})
-  x, iterator = x_and_iter
-  it, x_star = ProximalAlgorithms.run!(iterator)
-  ~x .= x_star
-  return it, iterator 
+  error("Sorry, I cannot parse this problem for solver of type $(T)")
 end
 
 
 export solve
 
 """
-`solve(terms::Tuple, solver_opt::ForwardBackwardSolver)`
+    solve(terms::Tuple, solver::ForwardBackwardSolver)
 
 Takes as input a tuple containing the terms defining the problem and the solver options.
 
@@ -102,22 +58,26 @@ Solves the problem returning a tuple containing the iterations taken and the bui
 # Example
 
 ```julia
-
 julia> x = Variable(4)
 Variable(Float64, (4,))
 
 julia> A, b = randn(10,4), randn(10);
 
-julia> solve(p,PG());
-it |      gamma |        fpr |
-------|------------|------------|
-1 | 7.6375e-02 | 1.8690e+00 |
-12 | 7.6375e-02 | 9.7599e-05 |
+julia> p = problem(ls(A*x - b ), norm(x) <= 1);
 
-```
+julia> solve(p, ProximalAlgorithms.ForwardBackward());
 
+julia> ~x
+4-element Array{Float64,1}:
+ -0.6427139974173074
+ -0.29043653211431103
+ -0.6090539651510192
+  0.36279278640995494
+```
 """
 function solve(terms::Tuple, solver::ForwardBackwardSolver)
-  built_slv = build(terms, solver)
-  return solve!(built_slv)
+    x, kwargs = parse_problem(terms, solver)
+    x_star, it = solver(~x; kwargs...)
+    ~x .= x_star
+    return x, it
 end

src/solvers/solvers_options.jl

@@ -1,105 +1,9 @@
-abstract type Solver end
+using ProximalAlgorithms
 
-abstract type ForwardBackwardSolver <: Solver end
+const ForwardBackwardSolver = Union{
+    ProximalAlgorithms.ForwardBackward,
+    ProximalAlgorithms.ZeroFPR,
+    ProximalAlgorithms.PANOC,
+}
 
-export Solver, ForwardBackwardSolver
-
-################################################################################
-export PG, FPG
-
-"""
-`PG(;kwargs...)`
-
-Creates an object `PG` containing the options of the Proximal Gradient solvers:
-
-  * `gamma`, stepsize (default: unspecified, determined automatically)
-  * `maxit`, maximum number of iteration (default: `10000`)
-  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
-  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
-  * `fast`, enables accelerated method (default: `false`)
-  * `verbose`, verbosity level (default: `1`)
-  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)
-
-"""
-struct PG <: ForwardBackwardSolver
-    kwargs::Iterators.Pairs
-    function PG(; kwargs...)
-        new(kwargs)
-    end
-end
-
-"""
-`FPG(;kwargs...)`
-
-Same as `PG`, creates the options of the Fast Proximal Gradient solver. 
-
-"""
-function FPG(; kwargs...)
-    return PG(; kwargs..., fast=true)
-end
-
-function build_iterator(x, solver::PG; kwargs...)
-    x, ProximalAlgorithms.FBSIterator(~x; solver.kwargs..., kwargs...)
-end
-
-################################################################################
-export ZeroFPR
-
-"""
-`ZeroFPR(;kwargs...)`
-
-Creates an object `ZeroFPR` containing the options of the ZeroFPR solver:
-
-  * `gamma`, stepsize (default: unspecified, determined automatically)
-  * `maxit`, maximum number of iteration (default: `10000`)
-  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
-  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
-  * `fast`, enables accelerated method (default: `false`)
-  * `verbose`, verbosity level (default: `1`)
-  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)
-  * `memory`, memory of the `LBFGS` operator (default: `10` )
-
-"""
-struct ZeroFPR <: ForwardBackwardSolver
-    kwargs::Iterators.Pairs
-    function ZeroFPR(; kwargs...)
-        new(kwargs)
-    end
-end
-
-function build_iterator(x, solver::ZeroFPR; kwargs...)
-    x, ProximalAlgorithms.ZeroFPRIterator(~x; solver.kwargs..., kwargs...)
-end
-
-################################################################################
-export PANOC
-
-"""
-`ZeroFPR(;kwargs...)`
-
-Creates an object `PANOC` containing the options of the PANOC solver:
-
-  * `gamma`, stepsize (default: unspecified, determined automatically)
-  * `maxit`, maximum number of iteration (default: `10000`)
-  * `tol`, halting tolerance on the fixed-point residual (default: `1e-4`)
-  * `adaptive`, adaptively adjust `gamma` (default: `false` if `gamma` is provided)
-  * `fast`, enables accelerated method (default: `false`)
-  * `verbose`, verbosity level (default: `1`)
-  * `verbose_freq`, verbosity frequency for `verbose = 1` (default: `100`)
-  * `memory`, memory of the `LBFGS` operator (default: `10` )
-
-"""
-struct PANOC <: ForwardBackwardSolver
-    kwargs::Iterators.Pairs
-    function PANOC(; kwargs...)
-        new(kwargs)
-    end
-end
-
-function build_iterator(x, solver::PANOC; kwargs...)
-    x, ProximalAlgorithms.PANOCIterator(~x; solver.kwargs..., kwargs...)
-end
-
-default_solver = PANOC
-
-################################################################################
+const default_solver = ProximalAlgorithms.PANOC