StructuredOptimization-related changes + new algorithms

.github/workflows/ci.yml

@@ -6,28 +6,31 @@ on:
       - master
   pull_request:
   workflow_dispatch:
+
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }}
+    name: Julia ${{ matrix.julia-version }} - ${{ matrix.os }} - ${{ matrix.julia-arch }}
     runs-on: ${{ matrix.os }}
+
     strategy:
       fail-fast: false
       matrix:
-        version:
-          - '1'
-          - '1.6'
-        os:
-          - ubuntu-latest
-          - macOS-latest
-          - windows-latest
-        arch:
-          - x64
+        julia-version: ['1.6', 'lts', '1']
+        julia-arch: [x64]
+        os: [ubuntu-latest, windows-latest, macOS-latest]
+
+    # needed to allow julia-actions/cache to delete old caches that it has created
+    permissions:
+      actions: write
+      contents: read
+
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v5
       - uses: julia-actions/setup-julia@v1
         with:
-          version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
+          version: ${{ matrix.julia-version }}
+          arch: ${{ matrix.julia-arch }}
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
 

Project.toml

@@ -1,18 +1,43 @@
 name = "ProximalAlgorithms"
 uuid = "140ffc9f-1907-541a-a177-7475e0a401e9"
-version = "0.7.0"
+version = "0.8.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OperatorCore = "3945cd23-d97e-4db0-9df2-35342dbd287d"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProximalCore = "dc4f5ac2-75d1-4f31-931e-60435d74994b"
 
 [compat]
 ADTypes = "1.5.3"
-DifferentiationInterface = "0.6.2"
+AbstractOperators = "0.4"
+Aqua = "0.8"
+DifferentiationInterface = "0.6.2,0.7"
+ForwardDiff = "0.10"
 LinearAlgebra = "1.2"
+OperatorCore = "0.1"
 Printf = "1.2"
-ProximalCore = "0.1"
+ProximalCore = "0.2"
+ProximalOperators = "0.17"
+Random = "1"
+RecursiveArrayTools = "3.31"
+ReverseDiff = "1.15"
+Test = "1.11"
+Zygote = "0.7"
 julia = "1.6"
+
+[extras]
+AbstractOperators = "d9c5613a-d543-52d8-9afd-8f241a8c3f1c"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+
+[targets]
+test = ["AbstractOperators", "Aqua", "ForwardDiff", "ProximalOperators", "Random", "RecursiveArrayTools", "ReverseDiff", "Test", "Zygote"]

src/ProximalAlgorithms.jl

@@ -3,7 +3,15 @@ module ProximalAlgorithms
 using ADTypes: ADTypes
 using DifferentiationInterface: DifferentiationInterface
 using ProximalCore
-using ProximalCore: prox, prox!
+using ProximalCore: Zero, IndZero, convex_conjugate, prox, prox!, is_smooth, is_locally_smooth, is_convex, is_strongly_convex, is_proximable
+using OperatorCore: is_linear
+using LinearAlgebra
+using Base.Iterators
+using Printf
+
+import Base: show
+import Base: *
+import LinearAlgebra: mul!
 
 const RealOrComplex{R} = Union{R,Complex{R}}
 const Maybe{T} = Union{T,Nothing}
@@ -39,6 +47,22 @@ function value_and_gradient(f::ProximalCore.Zero, x)
     return f(x), zero(x)
 end
 
+"""
+    value_and_gradient!(grad_f_x, f, x)
+
+Compute the value of `f` at `x` and store the gradient in `grad_f_x`.
+Returns the value of `f` at `x`.
+"""
+function value_and_gradient!(grad_f_x, f::AutoDifferentiable, x)
+    f_x, grad_f_x = DifferentiationInterface.value_and_gradient!(f.f, grad_f_x, f.backend, x)
+    return f_x
+end
+
+function value_and_gradient!(grad_f_x, f::ProximalCore.Zero, x)
+    fill!(grad_f_x, 0)
+    return f(x)
+end
+
 # various utilities
 
 include("utilities/fb_tools.jl")
@@ -55,18 +79,19 @@ include("accel/noaccel.jl")
 
 # algorithm interface
 
-struct IterativeAlgorithm{IteratorType,H,S,D,K}
+struct IterativeAlgorithm{IteratorType,H,S,I,D,K}
     maxit::Int
     stop::H
     solution::S
     verbose::Bool
     freq::Int
+    summary::I
     display::D
     kwargs::K
 end
 
 """
-    IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...)
+    IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, summary, display, kwargs...)
 
 Wrapper for an iterator type `T`, adding termination and verbosity options on top of it.
 
@@ -75,7 +100,7 @@ The resulting "algorithm" object `alg` can be called on a set of keyword argumen
 to `kwargs` and passed on to `T` to construct an iterator which will be looped over.
 Specifically, if an algorithm is constructed as
 
-    alg = IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...)
+    alg = IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, summary, display, kwargs...)
 
 then calling it with
 
@@ -88,7 +113,7 @@ will internally loop over an iterator constructed as
 # Note
 This constructor is not meant to be used directly: instead, algorithm-specific constructors
 should be defined on top of it and exposed to the user, that set appropriate default functions
-for `stop`, `solution`, `display`.
+for `stop`, `solution`, `summary`, `display`.
 
 # Arguments
 * `T::Type`: iterator type to use
@@ -97,35 +122,144 @@ for `stop`, `solution`, `display`.
 * `solution::Function`: solution mapping, `solution(::T, state)` should return the identified solution
 * `verbose::Bool`: whether the algorithm state should be displayed
 * `freq::Int`: every how many iterations to display the algorithm state
-* `display::Function`: display function, `display(::Int, ::T, state)` should display a summary of the iteration state
+* `summary::Function`: function returning a summary of the iteration state, `summary(k::Int, iter::T, state)` should return a vector of pairs `(name, value)`
+* `display::Function`: display function, `display(k::Int, alg, iter::T, state)` should display a summary of the iteration state
 * `kwargs...`: keyword arguments to pass on to `T` when constructing the iterator
 """
-IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, display, kwargs...) =
-    IterativeAlgorithm{T,typeof(stop),typeof(solution),typeof(display),typeof(kwargs)}(
+IterativeAlgorithm(T; maxit, stop, solution, verbose, freq, summary, display, kwargs...) =
+    IterativeAlgorithm{T,typeof(stop),typeof(solution),typeof(summary),typeof(display),typeof(kwargs)}(
         maxit,
         stop,
         solution,
         verbose,
         freq,
+        summary,
         display,
         kwargs,
     )
 
+"""
+    override_parameters(alg::IterativeAlgorithm; new_kwargs...)
+
+Return a new `IterativeAlgorithm` of the same type as `alg`, but with parameters overridden by `new_kwargs`.
+This is a convenience function to allow for easy modification of an existing algorithm object.
+"""
+function override_parameters(alg::IterativeAlgorithm; new_kwargs...)
+    if isempty(new_kwargs)
+        return alg
+    end
+    kwargs = Dict{Symbol, Any}(
+        :maxit => alg.maxit,
+        :stop => alg.stop,
+        :solution => alg.solution,
+        :verbose => alg.verbose,
+        :freq => alg.freq,
+        :summary => alg.summary,
+        :display => alg.display)
+    merge!(kwargs, alg.kwargs)
+    merge!(kwargs, new_kwargs)
+    return IterativeAlgorithm(typeof(alg).parameters[1]; kwargs...)
+end
+
+"""
+    get_iterator(alg::IterativeAlgorithm{IteratorType}) where {IteratorType}
+
+Return an iterator of type `IteratorType` constructed from the algorithm `alg`.
+This is a convenience function to allow for easy access to the iterator type
+associated with an `IterativeAlgorithm`.
+
+# Example
+```julia
+julia> using ProximalAlgorithms: CG, get_iterator
+
+julia> alg = CG(maxit=3, tol=1e-8);
+
+julia> iter = get_iterator(alg, A=reshape(collect(1:25)), b=collect(1:5));
+
+julia>  for (k, state) in enumerate(iter)
+            if k >= alg.maxit || alg.stop(iter, state)
+                alg.verbose && alg.display(k, iter, state)
+                return (alg.solution(iter, state), k)
+            end
+            alg.verbose && mod(k, alg.freq) == 0 && alg.display(k, iter, state)
+        end
+    1 | 7.416e+00
+    2 | 2.742e+00
+    3 | 2.300e+01
+([0.5581699346405239, 0.31633986928104635, 0.07450980392156867, -0.16732026143790907, -0.4091503267973867], 3)
+```
+"""
+get_iterator(alg::IterativeAlgorithm{IteratorType}; kwargs...) where {IteratorType} =
+    IteratorType(; alg.kwargs..., kwargs...)
+
+function default_display(k, alg, iter, state, printfunc=println)
+    if alg.freq > 0
+        summary = alg.summary(k, iter, state)
+        column_widths = map(pair -> max(length(pair.first), pair.second isa Integer ? 5 : 9), summary)
+        if k == 0
+            keys = map(first, summary)
+            first_line = [_get_centered_text(key, width) for (width, key) in zip(column_widths, keys)]
+            printfunc(join(first_line, " | "))
+            second_line = [repeat('-', width) for width in column_widths]
+            printfunc(join(second_line, "-|-"), "-")
+        else
+            values = map(last, summary)
+            parts = [_format_value(value, width) for (width, value) in zip(column_widths, values)]
+            printfunc(join(parts, " | "))
+        end
+    else
+        summary = alg.summary(k, iter, state)
+        if summary[1].first == ""
+            summary = ("total iterations" => k, summary[2:end]...)
+        end
+        items = map(pair -> @sprintf("%s=%s", pair.first, _format_value(pair.second, 0)), summary)
+        printfunc(join(items, ", "))
+    end
+end
+
+function _get_centered_text(text, width)
+    l = length(text)
+    if l >= width
+        return text
+    end
+    left_padding = div(width - l, 2)
+    right_padding = width - l - left_padding
+    return repeat(" ", left_padding) * text * repeat(" ", right_padding)
+end
+
+function _format_value(value, width)
+    if value isa Integer
+        return @sprintf("%*d", width, value)
+    elseif value isa Float64 || value isa Float32
+        return @sprintf("%*.3e", width, value)
+    else
+        return @sprintf("%*s", width, string(value))
+    end
+end
+
 function (alg::IterativeAlgorithm{IteratorType})(; kwargs...) where {IteratorType}
     iter = IteratorType(; alg.kwargs..., kwargs...)
     for (k, state) in enumerate(iter)
+        if k == 1 && alg.verbose && alg.freq > 0
+            alg.display(0, alg, iter, state)
+        end
         if k >= alg.maxit || alg.stop(iter, state)
-            alg.verbose && alg.display(k, iter, state)
+            alg.verbose && alg.display(k, alg, iter, state)
             return (alg.solution(iter, state), k)
         end
-        alg.verbose && mod(k, alg.freq) == 0 && alg.display(k, iter, state)
+        alg.verbose && alg.freq > 0 && mod(k, alg.freq) == 0 && alg.display(k, alg, iter, state)
     end
 end
 
+include("utilities/get_assumptions.jl")
+
 # algorithm implementations
 
+include("algorithms/cg.jl")
+include("algorithms/admm.jl")
 include("algorithms/forward_backward.jl")
 include("algorithms/fast_forward_backward.jl")
+include("algorithms/pogm.jl")
 include("algorithms/zerofpr.jl")
 include("algorithms/panoc.jl")
 include("algorithms/douglas_rachford.jl")
@@ -136,4 +270,31 @@ include("algorithms/li_lin.jl")
 include("algorithms/sfista.jl")
 include("algorithms/panocplus.jl")
 
+include("penalty_sequences/penalty_sequence_base.jl")
+include("penalty_sequences/fixed_penalty.jl")
+include("penalty_sequences/residual_balancing_penalty.jl")
+include("penalty_sequences/wohlberg_penalty.jl")
+include("penalty_sequences/barzilai_borwein_penalty.jl")
+include("penalty_sequences/spectral_radius_approx_penalty.jl")
+include("penalty_sequences/spectral_radius_bound_penalty.jl")
+
+get_algorithms() = [
+    CG(),
+    CGNR(),
+    FastForwardBackward(),
+    POGM(),
+    ZeroFPR(),
+    PANOCplus(),
+    DavisYin(),
+    VuCondat(),
+    DouglasRachford(),
+    ADMM(),
+    SFISTA(),
+    DRLS(),
+    ChambollePock(),
+    LiLin(),
+    PANOC(),
+    ForwardBackward(),
+]
+
 end # module

src/accel/anderson.jl

@@ -1,7 +1,3 @@
-using LinearAlgebra
-import Base: *
-import LinearAlgebra: mul!
-
 mutable struct AndersonAccelerationOperator{M,I,T}
     currmem::I
     curridx::I

src/accel/broyden.jl

@@ -1,7 +1,3 @@
-using LinearAlgebra
-import Base: *
-import LinearAlgebra: mul!
-
 struct BroydenOperator{R,TH}
     H::TH
     theta_bar::R

src/accel/lbfgs.jl

@@ -1,7 +1,3 @@
-using LinearAlgebra
-import Base: *
-import LinearAlgebra: mul!
-
 mutable struct LBFGSOperator{M,R,I,T}
     currmem::I
     curridx::I