Proximal Point Method (#410)

* Update the ProximalProblem to also cover single function proxes.
* Implement a proximal point state.
* Implement and document PPA.
* bump version, changelog.
* improve get proximal map – add test coverage.
* Finish test coverage and changelog.
* a vale check.
* reintroduce code cov for counting single prox access.

Changelog.md

@@ -5,6 +5,16 @@ All notable Changes to the Julia package `Manopt.jl` will be documented in this
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+# [0.5.1] – September 4, 2024
+
+## Changed
+
+* slightly improves the test for the ` ExponentialFamilyProjection` text on the about page.
+
+## Added
+
+* the `proximal_point` method.
+
 # [0.5.0] – August 29, 2024
 
 This breaking update is mainly concerned with improving a unified experience through all solvers

Project.toml

@@ -1,7 +1,7 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 authors = ["Ronny Bergmann <[email protected]>"]
-version = "0.5.0"
+version = "0.5.1"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"

docs/src/about.md

@@ -12,7 +12,7 @@ Thanks to the following contributors to `Manopt.jl`:
 * [Willem Diepeveen](https://www.maths.cam.ac.uk/person/wd292) implemented the [primal-dual Riemannian semismooth Newton](solvers/primal_dual_semismooth_Newton.md) solver.
 * [Hajg Jasa](https://www.ntnu.edu/employees/hajg.jasa) implemented the [convex bundle method](solvers/convex_bundle_method.md) and the [proximal bundle method](solvers/proximal_bundle_method.md) and a default subsolver each of them.
 * Even Stephansen Kjemsås contributed to the implementation of the [Frank Wolfe Method](solvers/FrankWolfe.md) solver.
-* Mathias Ravn Munkvold contributed most of the implementation of the [Adaptive Regularization with Cubics](solvers/adaptive-regularization-with-cubics.md) solver as well as ist [Lanczos](@ref arc-Lanczos) subsolver
+* Mathias Ravn Munkvold contributed most of the implementation of the [Adaptive Regularization with Cubics](solvers/adaptive-regularization-with-cubics.md) solver as well as its [Lanczos](@ref arc-Lanczos) subsolver
 * [Tom-Christian Riemer](https://www.tu-chemnitz.de/mathematik/wire/mitarbeiter.php) implemented the [trust regions](solvers/trust_regions.md) and [quasi Newton](solvers/quasi_Newton.md) solvers as well as the [truncated conjugate gradient descent](solvers/truncated_conjugate_gradient_descent.md) subsolver.
 * [Markus A. Stokkenes](https://www.linkedin.com/in/markus-a-stokkenes-b41bba17b/) contributed most of the implementation of the [Interior Point Newton Method](solvers/interior_point_Newton.md) as well as its default [Conjugate Residual](solvers/conjugate_residual.md) subsolver
 * [Manuel Weiss](https://scoop.iwr.uni-heidelberg.de/author/manuel-weiß/) implemented most of the [conjugate gradient update rules](@ref cg-coeffs)

docs/src/plans/stepsize.md

@@ -37,7 +37,7 @@ Tangent bundle with the Sasaki metric has 0 injectivity radius, so the maximum s
 `Hyperrectangle` also has 0 injectivity radius and an estimate based on maximum of dimensions along each index is used instead.
 For manifolds with corners, however, a line search capable of handling break points along the projected search direction should be used, and such algorithms do not call `max_stepsize`.
 
-Internally the step size functions above create a [`ManifoldDefaultsFactory`](@ref).
+Internally these step size functions create a [`ManifoldDefaultsFactory`](@ref).
 Internally these use
 
 ```@autodocs

docs/src/solvers/conjugate_gradient_descent.md

@@ -32,7 +32,7 @@ PolakRibiereCoefficient
 SteepestDescentCoefficient
 ```
 
-## Internal rule storages
+## Internal rules for coefficients
 
 ```@docs
 Manopt.ConjugateGradientBealeRestartRule

docs/src/solvers/index.md

@@ -72,6 +72,7 @@ If the gradient does not exist everywhere, that is if the splitting yields summa
 * [Difference of Convex Proximal Point](@ref solver-difference-of-convex-proximal-point) uses a splitting of the (non-convex) function ``f = g - h`` into a difference of two functions; provided the proximal map of ``g`` and the subgradient of ``h``, the next iterate is computed. Compared to DCA, the corresponding sub problem is here written in a form that yields the proximal map.
 * [Douglas—Rachford](DouglasRachford.md) uses a splitting ``f(p) = F(x) + G(x)`` and their proximal maps to compute a minimizer of ``f``, which can be non-smooth.
 * [Primal-dual Riemannian semismooth Newton Algorithm](@ref solver-pdrssn) extends Chambolle-Pock and requires the differentials of the proximal maps additionally.
+* The [Proximal Point](proximal_point.md) uses the proximal map of ``f`` iteratively.
 
 ## Constrained
 
@@ -120,6 +121,7 @@ For these you can use
 | [Particle Swarm](particle_swarm.md) | [`particle_swarm`](@ref) | [`ParticleSwarmState`](@ref) |
 [Primal-dual Riemannian semismooth Newton Algorithm](@ref solver-pdrssn) | [`primal_dual_semismooth_Newton`](@ref) | [`PrimalDualSemismoothNewtonState`](@ref) |
 | [Proximal Bundle Method](proximal_bundle_method.md) | [`proximal_bundle_method`](@ref) | [`ProximalBundleMethodState`](@ref) |
+| [Proximal Point](proximal_point.md) | [`proximal_point`](@ref) |  [`ProximalPointState`](@ref) |
 | [Quasi-Newton Method](quasi_Newton.md) | [`quasi_Newton`](@ref) | [`QuasiNewtonState`](@ref) |
 | [Steihaug-Toint Truncated Conjugate-Gradient Method](@ref tCG) | [`truncated_conjugate_gradient_descent`](@ref) | [`TruncatedConjugateGradientState`](@ref) |
 | [Subgradient Method](subgradient.md) | [`subgradient_method`](@ref) | [`SubGradientMethodState`](@ref) |

docs/src/solvers/proximal_point.md

@@ -0,0 +1,21 @@
+# Proximal point method
+
+```@meta
+CurrentModule = Manopt
+```
+
+```@docs
+proximal_point
+proximal_point!
+```
+
+## State
+
+```@docs
+ProximalPointState
+```
+
+```@bibliography
+Pages = ["proximal_point.md"]
+Canonical=false
+```

src/Manopt.jl

@@ -201,6 +201,7 @@ include("solvers/LevenbergMarquardt.jl")
 include("solvers/particle_swarm.jl")
 include("solvers/primal_dual_semismooth_Newton.jl")
 include("solvers/proximal_bundle_method.jl")
+include("solvers/proximal_point.jl")
 include("solvers/quasi_Newton.jl")
 include("solvers/truncated_conjugate_gradient_descent.jl")
 include("solvers/trust_regions.jl")
@@ -507,6 +508,8 @@ export adaptive_regularization_with_cubics,
     primal_dual_semismooth_Newton,
     proximal_bundle_method,
     proximal_bundle_method!,
+    proximal_point,
+    proximal_point!,
     quasi_Newton,
     quasi_Newton!,
     stochastic_gradient_descent,

src/plans/proximal_plan.jl

@@ -6,24 +6,31 @@
 @doc """
     ManifoldProximalMapObjective{E<:AbstractEvaluationType, TC, TP, V <: Vector{<:Integer}} <: AbstractManifoldCostObjective{E, TC}
 
-specify a problem for solvers based on the evaluation of proximal maps.
+specify a problem for solvers based on the evaluation of proximal maps,
+which represents proximal maps ``$(_tex(:prox))_{λf_i}`` for summands ``f = f_1 + f_2+ … + f_N`` of the cost function ``f``.
 
 # Fields
 
-* `cost`: a function ``F:$(_tex(:Cal, "M"))→ℝ`` to
+* `cost`: a function ``f:$(_tex(:Cal, "M"))→ℝ`` to
   minimize
-* `proxes`: proximal maps ``$(_tex(:prox))_{λφ}:$(_tex(:Cal, "M")) → $(_tex(:Cal, "M"))``
-  as functions `(M, λ, p) -> q`.
+* `proxes`: proximal maps ``$(_tex(:prox))_{λf_i}:$(_tex(:Cal, "M")) → $(_tex(:Cal, "M"))``
+  as functions `(M, λ, p) -> q` or in-place `(M, q, λ, p)`.
 * `number_of_proxes`: number of proximal maps per function,
   to specify when one of the maps is a combined one such that the proximal maps
   functions return more than one entry per function, you have to adapt this value.
   if not specified, it is set to one prox per function.
 
 # Constructor
 
-    ManifoldProximalMapObjective(cost, proxes, numer_of_proxes=onex(length(proxes));
+    ManifoldProximalMapObjective(f, proxes_f::Union{Tuple,AbstractVector}, numer_of_proxes=onex(length(proxes));
        evaluation=Allocating)
 
+Generate a proximal problem with a tuple or vector of funtions, where by default every function computes a single prox
+of one component of ``f``.
+
+    ManifoldProximalMapObjective(f, prox_f); evaluation=Allocating)
+
+Generate a proximal objective for ``f`` and its proxial map ``$(_tex(:prox))_{λf}``
 
 # See also
 
@@ -60,10 +67,12 @@ mutable struct ManifoldProximalMapObjective{E<:AbstractEvaluationType,TC,TP,V} <
             new{E,F,typeof(proxes_f),typeof(nOP)}(f, proxes_f, nOP)
         end
     end
-end
-function check_prox_number(n, i)
-    (i > n) && throw(ErrorException("the $(i)th entry does not exists, only $n available."))
-    return true
+    function ManifoldProximalMapObjective(
+        f::F, prox_f::PF; evaluation::E=AllocatingEvaluation()
+    ) where {E<:AbstractEvaluationType,F,PF}
+        i = 1
+        return new{E,F,PF,typeof(i)}(f, prox_f, i)
+    end
 end
 @doc raw"""
     q = get_proximal_map(M::AbstractManifold, mpo::ManifoldProximalMapObjective, λ, p)
@@ -88,47 +97,97 @@ function get_proximal_map!(amp::AbstractManoptProblem, q, λ, p)
     return get_proximal_map!(get_manifold(amp), q, get_objective(amp), λ, p)
 end
 
+function check_prox_number(pf::Union{Tuple,Vector}, i)
+    n = length(pf)
+    (i > n) && throw(ErrorException("the $(i)th entry does not exists, only $n available."))
+    return true
+end
+
 function get_proximal_map(
-    M::AbstractManifold, mpo::ManifoldProximalMapObjective{AllocatingEvaluation}, λ, p, i
-)
-    check_prox_number(length(mpo.proximal_maps!!), i)
+    M::AbstractManifold,
+    mpo::ManifoldProximalMapObjective{AllocatingEvaluation,F,<:Union{<:Tuple,<:Vector}},
+    λ,
+    p,
+    i,
+) where {F}
+    check_prox_number(mpo.proximal_maps!!, i)
     return mpo.proximal_maps!![i](M, λ, p)
 end
 function get_proximal_map(
-    M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, λ, p, i
+    M::AbstractManifold, admo::AbstractDecoratedManifoldObjective, args...
 )
-    return get_proximal_map(M, get_objective(admo, false), λ, p, i)
+    return get_proximal_map(M, get_objective(admo, false), args...)
 end
-
 function get_proximal_map!(
-    M::AbstractManifold, q, mpo::ManifoldProximalMapObjective{AllocatingEvaluation}, λ, p, i
-)
-    check_prox_number(length(mpo.proximal_maps!!), i)
+    M::AbstractManifold,
+    q,
+    mpo::ManifoldProximalMapObjective{AllocatingEvaluation,F,<:Union{<:Tuple,<:Vector}},
+    λ,
+    p,
+    i,
+) where {F}
+    check_prox_number(mpo.proximal_maps!!, i)
     copyto!(M, q, mpo.proximal_maps!![i](M, λ, p))
     return q
 end
 function get_proximal_map!(
-    M::AbstractManifold, q, admo::AbstractDecoratedManifoldObjective, λ, p, i
+    M::AbstractManifold, q, admo::AbstractDecoratedManifoldObjective, args...
 )
-    return get_proximal_map!(M, q, get_objective(admo, false), λ, p, i)
+    return get_proximal_map!(M, q, get_objective(admo, false), args...)
 end
 function get_proximal_map(
-    M::AbstractManifold, mpo::ManifoldProximalMapObjective{InplaceEvaluation}, λ, p, i
-)
-    check_prox_number(length(mpo.proximal_maps!!), i)
+    M::AbstractManifold,
+    mpo::ManifoldProximalMapObjective{InplaceEvaluation,F,<:Union{<:Tuple,<:Vector}},
+    λ,
+    p,
+    i,
+) where {F}
+    check_prox_number(mpo.proximal_maps!!, i)
     q = allocate_result(M, get_proximal_map, p)
     mpo.proximal_maps!![i](M, q, λ, p)
     return q
 end
 function get_proximal_map!(
-    M::AbstractManifold, q, mpo::ManifoldProximalMapObjective{InplaceEvaluation}, λ, p, i
-)
-    check_prox_number(length(mpo.proximal_maps!!), i)
+    M::AbstractManifold,
+    q,
+    mpo::ManifoldProximalMapObjective{InplaceEvaluation,F,<:Union{<:Tuple,<:Vector}},
+    λ,
+    p,
+    i,
+) where {F}
+    check_prox_number(mpo.proximal_maps!!, i)
     mpo.proximal_maps!![i](M, q, λ, p)
     return q
 end
 #
 #
+# Single function accessors
+function get_proximal_map(
+    M::AbstractManifold, mpo::ManifoldProximalMapObjective{AllocatingEvaluation}, λ, p
+)
+    return mpo.proximal_maps!!(M, λ, p)
+end
+function get_proximal_map!(
+    M::AbstractManifold, q, mpo::ManifoldProximalMapObjective{AllocatingEvaluation}, λ, p
+)
+    copyto!(M, q, mpo.proximal_maps!!(M, λ, p))
+    return q
+end
+function get_proximal_map(
+    M::AbstractManifold, mpo::ManifoldProximalMapObjective{InplaceEvaluation}, λ, p
+)
+    q = allocate_result(M, get_proximal_map, p)
+    mpo.proximal_maps!!(M, q, λ, p)
+    return q
+end
+function get_proximal_map!(
+    M::AbstractManifold, q, mpo::ManifoldProximalMapObjective{InplaceEvaluation}, λ, p
+)
+    mpo.proximal_maps!!(M, q, λ, p)
+    return q
+end
+#
+#
 # Proximal based State
 #
 #
@@ -141,16 +200,20 @@ stores options for the [`cyclic_proximal_point`](@ref) algorithm. These are the
 
 $(_var(:Field, :p; add=[:as_Iterate]))
 $(_var(:Field, :stopping_criterion, "stop"))
-* `λ`:         a function for the values of ``λ_k`` per iteration(cycle ``ì``
+* `λ`:         a function for the values of ``λ_k`` per iteration(cycle ``k``
 * `oder_type`: whether to use a randomly permuted sequence (`:FixedRandomOrder`),
   a per cycle permuted sequence (`:RandomOrder`) or the default linear one.
 
 # Constructor
 
-    CyclicProximalPointState(M; kwargs...)
+    CyclicProximalPointState(M::AbstractManifold; kwargs...)
 
 Generate the options
 
+## Input
+
+$(_var(:Argument, :M; type=true))
+
 # Keyword arguments
 
 * `evaluation_order=:LinearOrder`: soecify the `order_type`

src/plans/record.jl

@@ -523,7 +523,7 @@ $(_var(:Keyword, :inverse_retraction_method))
         storage                   = StoreStateAction(M; store_points=Tuple{:Iterate})
     )
 
-with the preceding fields as keywords. For the `DefaultManifold` only the field storage is used.
+with the previous fields as keywords. For the `DefaultManifold` only the field storage is used.
 Providing the actual manifold moves the default storage to the efficient point storage.
 """
 mutable struct RecordChange{

src/solvers/LevenbergMarquardt.jl

@@ -19,7 +19,7 @@ The second signature performs the optimization in-place of `p`.
 # Input
 
 $(_var(:Argument, :M; type=true))
-* `f`:              a cost function ``f: $(_math(:M)) M→ℝ^d``
+* `f`:              a cost function ``f: $(_math(:M))→ℝ^d``
 * `jacobian_f`:     the Jacobian of ``f``. The Jacobian is supposed to accept a keyword argument
   `basis_domain` which specifies basis of the tangent space at a given point in which the
   Jacobian is to be calculated. By default it should be the `DefaultOrthonormalBasis`.

src/solvers/cyclic_proximal_point.jl

@@ -26,7 +26,7 @@ perform a cyclic proximal point algorithm. This can be done in-place of `p`.
 # Input
 
 $(_var(:Argument, :M; type=true))
-* `f`:        a cost function ``f: $(_math(:M)) M→ℝ`` to minimize
+* `f`:        a cost function ``f: $(_math(:M))→ℝ`` to minimize
 * `proxes_f`: an Array of proximal maps (`Function`s) `(M,λ,p) -> q` or `(M, q, λ, p) -> q` for the summands of ``f`` (see `evaluation`)
 
 where `f` and the proximal maps `proxes_f` can also be given directly as a [`ManifoldProximalMapObjective`](@ref) `mpo`