Introduce proxes. (#37)

* Introduce proxes
* Unify docs, introduce Changelog, finish testing proxes (adapting tests formerly done in Manopt).
* Unify variables and docs.

---------

Co-authored-by: Mateusz Baran <[email protected]>

.github/workflows/changelog.yml

@@ -0,0 +1,12 @@
+name: Check Changelog
+on:
+  pull_request:
+
+jobs:
+  Check-Changelog:
+    name: Check Changelog Action
+    runs-on: ubuntu-latest
+    steps:
+      - uses: tarides/changelog-check-action@v2
+        with:
+          changelog: Changelog.md

Changelog.md

@@ -0,0 +1,14 @@
+# Changelog
+
+All notable changes to this Julia package will be documented in this file.
+
+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.3.9] - December 8, 2023
+
+### Added
+
+* proximal map of the distance function, `prox_distance`, and its mutating variant, `prox_distance!`
+* this changelog
+* A GitHub Action to check that the Changelog is updated on every PR.

Project.toml

@@ -1,7 +1,7 @@
 name = "ManifoldDiff"
 uuid = "af67fdf4-a580-4b9f-bbec-742ef357defd"
 authors = ["Seth Axen <[email protected]>", "Mateusz Baran <[email protected]>", "Ronny Bergmann <[email protected]>"]
-version = "0.3.8"
+version = "0.3.9"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/make.jl

@@ -1,3 +1,21 @@
+#!/usr/bin/env julia
+#
+#
+
+if "--help" ∈ ARGS
+    println(
+        """
+     docs/make.jl
+
+Render the `ManifoldDiff.jl` documentation with optional arguments
+
+Arguments
+* `--help`         - print this help and exit without rendering the documentation
+* `--prettyurls`   – toggle the prettyurls part to true (which is otherwise only true on CI)
+""",
+    )
+    exit(0)
+end
 # (a) if docs is not the current active environment, switch to it
 # (from https://github.com/JuliaIO/HDF5.jl/pull/1020/) 
 if Base.active_project() != joinpath(@__DIR__, "Project.toml")
@@ -8,15 +26,15 @@ if Base.active_project() != joinpath(@__DIR__, "Project.toml")
     Pkg.instantiate()
 end
 
-using ManifoldDiff, ManifoldsBase, ManifoldDiff
+using ManifoldsBase, ManifoldDiff
 using Documenter, DocumenterCitations
 using FiniteDiff, ForwardDiff, ReverseDiff, FiniteDifferences, Zygote
 
 bib = CitationBibliography(joinpath(@__DIR__, "src", "references.bib"); style = :alpha)
 makedocs(;
-    format = Documenter.HTML(
-        prettyurls = get(ENV, "CI", nothing) == "true",
-        assets = ["assets/favicon.ico"],
+    format = Documenter.HTML(;
+        prettyurls = (get(ENV, "CI", nothing) == "true") || ("--prettyurls" ∈ ARGS),
+        assets = ["assets/favicon.ico", "assets/citations.css"],
     ),
     modules = [
         ManifoldDiff,

docs/src/assets/citations.css

@@ -0,0 +1,19 @@
+/* Taken from https://juliadocs.org/DocumenterCitations.jl/v1.2/styling/ */
+
+.citation dl {
+  display: grid;
+  grid-template-columns: max-content auto; }
+.citation dt {
+  grid-column-start: 1; }
+.citation dd {
+  grid-column-start: 2;
+  margin-bottom: 0.75em; }
+.citation ul {
+ padding: 0 0 2.25em 0;
+ margin: 0;
+ list-style: none;}
+.citation ul li {
+ text-indent: -2.25em;
+ margin: 0.33em 0.5em 0.5em 2.25em;}
+.citation ol li {
+ padding-left:0.75em;}

docs/src/library.md

@@ -51,3 +51,22 @@ Modules = [ManifoldDiff]
 Pages = ["riemannian_diff.jl"]
 Order = [:type, :function, :constant]
 ```
+
+## Proximal Maps
+
+Given a convex, lower semi-continuous function ``f\colon \mathcal M \to \mathbb R``, its proximal map is defined
+for some ``λ>0`` as [Bacak:2014](@cite)
+
+```math
+\operatorname{prox}_{λf}(p) := \operatorname*{arg\,min}_{q\in\mathcal M} \frac{1}{2λ}d^2_{\mathcal M}(p,q) + f(q).
+```
+
+Another name for the proximal map is _resolvent_.
+Intuitively this means to minimize the function ``f`` while at the same timme “staying close”
+to the argument ``p``.
+
+```@autodocs
+Modules = [ManifoldDiff]
+Pages = ["proximal_maps.jl"]
+Order = [:type, :function, :constant]
+```

docs/src/references.bib

@@ -1,11 +1,23 @@
 @book{AbsilMahonySepulchre:2008,
     AUTHOR    = {Absil, P.-A. and Mahony, R. and Sepulchre, R.},
     DOI       = {10.1515/9781400830244},
-    NOTE      = {available online at [press.princeton.edu/chapters/absil/](http://press.princeton.edu/chapters/absil/)},
+    NOTE      = {available online at \href{http://press.princeton.edu/chapters/absil/}{press.princeton.edu/chapters/absil/}},
     PUBLISHER = {Princeton University Press},
     TITLE     = {Optimization Algorithms on Matrix Manifolds},
     YEAR      = {2008},
 }
+
+@article{Bacak:2014,
+    AUTHOR  = {Bačák, M.},
+    DOI     = {10.1137/140953393},
+    JOURNAL = {SIAM Journal on Optimization},
+    NUMBER  = {3},
+    PAGES   = {1542--1566},
+    TITLE   = {Computing medians and means in Hadamard spaces},
+    VOLUME  = {24},
+    YEAR    = {2014}
+}
+
 @book{Boumal:2023,
     TITLE     = {An Introduction to Optimization on Smooth Manifolds},
     AUTHOR    = {Boumal, Nicolas},
@@ -14,7 +26,7 @@ @book{Boumal:2023
     EDITION   = {First},
     PUBLISHER = {Cambridge University Press},
     DOI       = {10.1017/9781009166164},
-    URLDATE   = {2023-07-13},
+    NOTE      = {Homepage to the book: \href{https://www.nicolasboumal.net/book/index.html}{nicolasboumal.net/book/index.html}.},
     ABSTRACT  = {Optimization on Riemannian manifolds-the result of smooth geometry and optimization merging into one elegant modern framework-spans many areas of science and engineering, including machine learning, computer vision, signal processing, dynamical systems and scientific computing. This text introduces the differential geometry and Riemannian geometry concepts that will help students and researchers in applied mathematics, computer science and engineering gain a firm mathematical grounding to use these tools confidently in their research. Its charts-last approach will prove more intuitive from an optimizer's viewpoint, and all definitions and theorems are motivated to build time-tested optimization algorithms. Starting from first principles, the text goes on to cover current research on topics including worst-case complexity and geodesic convexity. Readers will appreciate the tricks of the trade for conducting research and for numerical implementations sprinkled throughout the book.},
     ISBN      = {978-1-00-916616-4}
 }
@@ -28,4 +40,15 @@ @article{Zimmermann:2020
     PAGES         = {A2593-A2619},
     YEAR          = {2020},
     DOI           = {10.1137/19M1282878},
+}
+
+@article{WeinmannDemaretStorath:2014,
+    AUTHOR  = {Weinmann, Andreas and Demaret, Laurent and Storath, Martin},
+    DOI     = {10.1137/130951075},
+    JOURNAL = {SIAM Journal on Imaging Sciences},
+    NUMBER  = {4},
+    PAGES   = {2226--2257},
+    TITLE   = {Total variation regularization for manifold-valued data},
+    VOLUME  = {7},
+    YEAR    = {2014}
 }

src/ManifoldDiff.jl

@@ -171,6 +171,7 @@ include("derivatives.jl")
 include("differentials.jl")
 include("gradients.jl")
 include("Jacobi_fields.jl")
+include("proximal_maps.jl")
 include("subgradients.jl")
 
 include("riemannian_diff.jl")

src/proximal_maps.jl

@@ -0,0 +1,50 @@
+@doc raw"""
+    y = prox_distance(M::AbstractManifold, λ::Real, p_data, p [, r=2])
+    prox_distance!(M::AbstractManifold, q, λ::Real, p_data, p [, r=2])
+
+Compute the proximal map ``\operatorname{prox}_{λf}`` with
+parameter λ of ``f(p) = \frac{1}{r}d_{\mathcal M}^r(p_{\mathrm{data}},p)``.
+For the in-place variant the computation is done in place of `q`.
+
+# Input
+* `M`      a manifold `M`
+* `λ`      the prox parameter, a positive real number.
+* `p_data` a point on `M`.
+* `p`      the argument of the proximal map
+* `r`      (`2`) exponent of the distance.
+
+# Output
+* `q` – the result of the proximal map of ``f``
+
+For more details see [WeinmannDemaretStorath:2014](@cite)
+"""
+function prox_distance(M::AbstractManifold, λ::Real, p_data, p, r::Int = 2)
+    d = distance(M, p_data, p)
+    if r == 2
+        t = λ / (1 + λ)
+    elseif r == 1
+        t = (λ < d) ? λ / d : 1.0
+    else
+        throw(
+            ErrorException(
+                "Proximal Map of distance(M, p_data, p) not implemented for an exponent $(r) (requires 1 or 2)",
+            ),
+        )
+    end
+    return shortest_geodesic(M, p, p_data, t)
+end
+function prox_distance!(M::AbstractManifold, q, λ::Real, p_data, p, r::Int = 2)
+    d = distance(M, p_data, p)
+    if r == 2
+        t = λ / (1 + λ)
+    elseif r == 1
+        t = (λ < d) ? λ / d : 1.0
+    else
+        throw(
+            ErrorException(
+                "Proximal Map of distance(M, p_data, p) not implemented for an exponent $(r) (requires 1 or 2)",
+            ),
+        )
+    end
+    return shortest_geodesic!(M, q, p, p_data, t)
+end

test/runtests.jl

@@ -15,5 +15,6 @@ using RecursiveArrayTools
     include("manifold_specializations.jl")
     include("test_gradients.jl")
     include("test_derivatives.jl")
+    include("test_proximal_maps.jl")
     include("test_subgradients.jl")
 end

test/test_proximal_maps.jl

@@ -0,0 +1,26 @@
+using Manifolds, ManifoldDiff, Test
+
+@testset "proximal maps" begin
+    #
+    # Distance
+    p = [1.0, 0.0, 0.0]
+    q = [0.0, 1.0, 0.0]
+    M = Sphere(2)
+    # Error for r != 1 or 2
+    @test_throws ErrorException ManifoldDiff.prox_distance(M, 1.0, p, q, 3)
+    @test_throws ErrorException ManifoldDiff.prox_distance!(M, p, 1.0, p, q, 3)
+    # r = 1
+    @test distance(
+        M,
+        ManifoldDiff.prox_distance(M, distance(M, p, q) / 2, p, q, 1),
+        shortest_geodesic(M, p, q, 0.5),
+    ) < eps()
+    t = similar(p)
+    ManifoldDiff.prox_distance!(M, t, distance(M, p, q) / 2, p, q, 1)
+    @test t == ManifoldDiff.prox_distance(M, distance(M, p, q) / 2, p, q, 1)
+    # r = 2
+    ManifoldDiff.prox_distance!(M, t, 1.0, p, q, 2)
+    @test t == ManifoldDiff.prox_distance(M, 1.0, p, q, 2)
+    @test t == ManifoldDiff.prox_distance(M, 1.0, p, q) # 2 is also the default
+    @test distance(M, t, shortest_geodesic(M, p, q, 0.5)) ≈ 0 atol = 1e-15
+end