Introduce the Riemannian Chambolle-Pock (#40)

* introduces a Riemannian Chambolle Poch plan.
* Introduces all examples from the paper
* generalises the original (l/e)RCPA to also run with retractions, inverse retractions and vector transport instead of just with exp/log and parallel transport

Author: kellertuer

.github/workflows/ci.yml

@@ -11,7 +11,7 @@ jobs:
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
-        julia-version: [1.0, 1.4]
+        julia-version: [1.0, 1.5]
         os: [ubuntu-latest, macOS-latest, windows-latest]
     steps:
       - uses: actions/checkout@v2
@@ -27,4 +27,4 @@ jobs:
           file: ./lcov.info
           name: codecov-umbrella
           fail_ci_if_error: false
-        if: ${{ matrix.julia-version == '1.4' && matrix.os =='ubuntu-latest' }}
+        if: ${{ matrix.julia-version == '1.5' && matrix.os =='ubuntu-latest' }}

.gitignore

@@ -4,15 +4,12 @@
 /docs/build/
 /docs/site/
 Tutorials/.ipynb_checkpoints
-src/examples/**/*.asy
-src/examples/**/*.png
-src/examples/**/*.csv
-src/examples/**/*.tex
-src/examples/**/*.pdf
-src/examples/**/*.mp4
-src/examples/**/*.jld2
+examples/**/*.asy
+examples/**/*.png
+examples/**/*.csv
+examples/**/*.tex
+examples/**/*.pdf
+examples/**/*.mp4
+examples/**/*.jld2
 docs/src/tutorials/*.md
 .vscode
-examples/Total_Variation/S2_TV
-examples/Total_Variation/Hn_TV
-examples/Total_Variation/S1_TV

Project.toml

@@ -1,7 +1,7 @@
 name = "Manopt"
 uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
 authors = ["Ronny Bergmann <manopt@ronnybergmann.net>"]
-version = "0.2.10"
+version = "0.2.11"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
@@ -22,7 +22,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 ColorSchemes = "3.5.0"
 ColorTypes = "0.9.1, 0.10"
 Colors = "0.11.2, 0.12"
-Manifolds = "0.4.1"
+Manifolds = "0.4.8"
 ManifoldsBase = "0.9"
 StaticArrays = "0.12"
 julia = "1.0"

docs/make.jl

@@ -18,13 +18,7 @@ for (i, tutorial) in enumerate(tutorials)
     global TutorialMenu
     sourceFile = joinpath(tutorialsInputPath, tutorial * ".jl")
     targetFile = joinpath(tutorialsOutputPath, tutorial * "md")
-    Literate.markdown(
-        sourceFile,
-        tutorialsOutputPath;
-        name=tutorial,
-        # codefence = "```julia" => "```",
-        credit=false,
-    )
+    Literate.markdown(sourceFile, tutorialsOutputPath; name=tutorial, credit=false)
     push!(TutorialMenu, menuEntries[i] => joinpath(tutorialsRelativePath, tutorial * ".md"))
 end
 makedocs(;
@@ -38,10 +32,11 @@ makedocs(;
         "Plans" => "plans/index.md",
         "Solvers" => [
             "Introduction" => "solvers/index.md",
+            "Chambolle-Pock" => "solvers/ChambollePock.md",
             "Conjugate gradient descent" => "solvers/conjugate_gradient_descent.md",
             "Cyclic Proximal Point" => "solvers/cyclic_proximal_point.md",
             "Douglas–Rachford" => "solvers/DouglasRachford.md",
-            "Gradient Descent" => "solvers/gradientDescent.md",
+            "Gradient Descent" => "solvers/gradient_descent.md",
             "Nelder–Mead" => "solvers/NelderMead.md",
             "Particle Swarm Optimization" => "solvers/particle_swarm.md",
             "Subgradient method" => "solvers/subgradient.md",
@@ -52,13 +47,13 @@ makedocs(;
         "Functions" => [
             "Introduction" => "functions/index.md",
             "Bézier curves" => "functions/bezier.md",
-            "Cost functions" => "functions/costFunctions.md",
+            "Cost functions" => "functions/costs.md",
             "Differentials" => "functions/differentials.md",
-            "Adjoint Differentials" => "functions/adjointDifferentials.md",
+            "Adjoint Differentials" => "functions/adjoint_differentials.md",
             "Gradients" => "functions/gradients.md",
-            "JacobiFields" => "functions/jacobiFields.md",
-            "Proximal Maps" => "functions/proximalMaps.md",
-            "Specific manifold functions" => "functions/manifold.md",
+            "Jacobi Fields" => "functions/Jacobi_fields.md",
+            "Proximal Maps" => "functions/proximal_maps.md",
+            "Specific Manifold Functions" => "functions/manifold.md",
         ],
         "Helpers" => [
             "Data" => "helpers/data.md",

docs/src/functions/jacobiFields.md renamed to docs/src/functions/Jacobi_fields.md

@@ -26,5 +26,5 @@ The following weights functions are available
 
 ```@autodocs
 Modules = [Manopt]
-Pages   = ["jacobiFields.jl"]
+Pages   = ["Jacobi_fields.jl"]
 ```

docs/src/functions/adjointDifferentials.md renamed to docs/src/functions/adjoint_differentials.md

@@ -2,5 +2,5 @@
 
 ```@autodocs
 Modules = [Manopt]
-Pages   = ["adjointDifferentials.jl"]
+Pages   = ["adjointdifferentials.jl"]
 ```

docs/src/functions/costFunctions.md renamed to docs/src/functions/costs.md

@@ -17,9 +17,9 @@ $\displaystyle\operatorname{prox}_{\lambda\varphi}(x^\star) = x^\star,$
 i.e. a minimizer is a fixed point of the proximal map.
 
 This page lists all proximal maps available within Manopt. To add you own, just
-extend the `functions/proximalMaps.jl` file.
+extend the `functions/proximal_maps.jl` file.
 
 ```@autodocs
 Modules = [Manopt]
-Pages   = ["proximalMaps.jl"]
+Pages   = ["proximal_maps.jl"]
 ```

docs/src/functions/proximalMaps.md renamed to docs/src/functions/proximal_maps.md

@@ -78,7 +78,7 @@ Further specific [`RecordAction`](@ref)s can be found at the specific Options.
 there's one internal helper that might be useful for you own actions, namely
 
 ```@docs
-record_or_eset!
+record_or_reset!
 ```
 
 ### [Stepsize and Linesearch](@id Stepsize)
@@ -141,7 +141,7 @@ get_subgradient
 
 ```@docs
 ProximalProblem
-getProximalMap
+get_proximal_map
 ```
 
 ### Further planned problems

docs/src/plans/index.md

@@ -0,0 +1,115 @@
+# [The Riemannian Chambolle-Pock Algorithm](@id ChambollePockSolver)
+
+The Riemannian Chambolle–Pock is a generalization of the Chambolle–Pock algorithm[^ChambollePock2011].
+It is also known as primal dual hybrig gradient (PDHG) or primal dual proximal splitting (PDPS) algorithm.
+
+In order to minimize over $p\in\mathcal M§ the cost function consisting of
+
+```math
+F(p) + G(\Lambda(p)),
+```
+
+where $F:\mathcal M \to \overline{\mathbb R}$, $G:\mathcal N \to \overline{\mathbb R}$, and
+$\Lambda:\mathcal M \to\mathcal N$.
+If the manifolds $\mathcal M$ or $\mathcal N$ are not Hadamard, it has to be considered locally,
+i.e. on geodesically convex sets $\mathcal C \subset \mathcal M$ and $\mathcal D \subset\mathcal N$
+such that $\Lambda(\mathcal C) \subset \mathcal D$.
+
+The algorithm is available in four variants: exact versus linearized (see `variant`)
+as well as with primal versus dual relaxation (see `relax`). For more details, see
+[^BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020].
+In the following we note the case of the exact, primal relaxed Riemannian Chambolle–Pock algorithm.
+
+Given base points $m\in\mathcal C$, $n=\Lambda(m)\in\mathcal D$,
+initial primal and dual values $p^{(0)} \in \mathcal C$, $\xi_n^{(0)} \in T_n^*\mathcal N$,
+and primal and dual step sizes $\sigma_0$, $\tau_0$, relaxation $\theta_0$,
+as well as acceleration $\gamma$.
+
+As an initialization, perform $\bar p^{(0)} \gets p^{(0)}$.
+
+The algorithms performs the steps $k=1,\ldots,$ (until a [`StoppingCriterion`](@ref) is fulfilled with)
+
+1. ```math
+   \xi^{(k+1)}_n = \operatorname{prox}_{\tau_k G_n^*}\Bigl(\xi_n^{(k)} + \tau_k \bigl(\log_n \Lambda (\bar p^{(k)})\bigr)^\flat\Bigr)
+   ```
+2. ```math
+   p^{(k+1)} = \operatorname{prox}_{\sigma_k F}\biggl(\exp_{p^{(k)}}\Bigl( \operatorname{PT}_{p^{(k)}\gets m}\bigl(-\sigma_k D\Lambda(m)^*[\xi_n^{(k+1)}]\bigr)^\sharp\Bigr)\biggr)
+   ```
+3. Update
+   * ``\theta_k = (1+2\gamma\sigma_k)^{-\frac{1}{2}}``
+   * ``\sigma_{k+1} = \sigma_k\theta_k``
+   * ``\tau_{k+1} =  \frac{\tau_k}{\theta_k}``
+4. ```math
+   \bar p^{(k+1)}  = \exp_{p^{(k+1)}}\bigl(-\theta_k \log_{p^{(k+1)}} p^{(k)}\bigr)
+   ```
+
+Furthermore you can exchange the exponential map, the logarithmic map, and the parallel transport
+by a retraction, an in verse retraction and a vector transport.
+
+Finally you can also update the base points $m$ and $n$ during the iterations.
+This introduces a few additional vector transports. The same holds for the case that
+$\Lambda(m^{(k)})\neq n^{(k)}$ at some point. All these cases are covered in the algorithm.
+
+```@meta
+CurrentModule = Manopt
+```
+
+```@docs
+ChambollePock
+```
+
+## Problem & Options
+
+```@docs
+PrimalDualProblem
+PrimalDualOptions
+ChambollePockOptions
+```
+
+## Useful Terms
+
+```@docs
+primal_residual
+dual_residual
+```
+
+## Debug
+
+```@docs
+DebugDualBaseIterate
+DebugDualBaseChange
+DebugPrimalBaseIterate
+DebugPrimalBaseChange
+DebugDualChange
+DebugDualIterate
+DebugDualResidual
+DebugPrimalChange
+DebugPrimalIterate
+DebugPrimalResidual
+DebugPrimalDualResidual
+```
+
+## Record
+
+```@docs
+RecordDualBaseIterate
+RecordDualBaseChange
+RecordDualChange
+RecordDualIterate
+RecordPrimalBaseIterate
+RecordPrimalBaseChange
+RecordPrimalChange
+RecordPrimalIterate
+```
+
+## Internals
+
+```@docs
+Manopt.update_prox_parameters!
+```
+
+[^ChambollePock2011]:
+    > A. Chambolle, T. Pock:
+    > _A first-order primal-dual algorithm for convex problems with applications to imaging_,
+    > Journal of Mathematical Imaging and Vision 40(1), 120–145, 2011.
+    > doi: [10.1007/s10851-010-0251-1](https://dx.doi.org/10.1007/s10851-010-0251-1)

docs/src/solvers/ChambollePock.md

@@ -14,10 +14,12 @@ The following algorithms are currently available
 
 | Solver  | File   | Problem & Option  |
 ----------|--------|-------------------|
-[steepest Descent](@ref GradientDescentSolver) | `gradient_descent.jl` |  [`GradientProblem`](@ref), [`GradientDescentOptions`](@ref)
 [Cyclic Proximal Point](@ref CPPSolver) | `cyclic_proximal_point.jl` | [`ProximalProblem`](@ref), [`CyclicProximalPointOptions`](@ref)
+[Chambolle-Pock](@ref ChambollePockSolver) | `Chambolle-Pock` | [`PrimalDualProblem`](@ref), [`ChambollePockOptions`](@ref)
 [Douglas–Rachford](@ref DRSolver) | `DouglasRachford.jl` | [`ProximalProblem`](@ref), [`DouglasRachfordOptions`](@ref)
+[Gradient Descent](@ref GradientDescentSolver) | `gradient_descent.jl` |  [`GradientProblem`](@ref), [`GradientDescentOptions`](@ref)
 [Nelder-Mead](@ref NelderMeadSolver) | `NelderMead.jl` | [`CostProblem`](@ref), [`NelderMeadOptions`](@ref)
+[Particle Swarm](@ref ParticleSwarmSolver) | `particle_swarm.jl` | [`CostProblem`](@ref), [`ParticleSwarmOptions`](@ref)
 [Subgradient Method](@ref SubgradientSolver) | `subgradient_method.jl` | [`SubGradientProblem`](@ref), [`SubGradientMethodOptions`](@ref)
 [Steihaug-Toint Truncated Conjugate-Gradient Method](@ref tCG) | `truncated_conjugate_gradient_descent.jl` | [`HessianProblem`](@ref),
 [`TruncatedConjugateGradientOptions`](@ref)

docs/src/solvers/gradientDescent.md renamed to docs/src/solvers/gradient_descent.md

@@ -27,7 +27,7 @@ $z_0 = \operatorname{P}(r_0)$, $\delta_0 = z_0$ and $k=0$
 
 Repeat until a convergence criterion is reached
 
-1. Set $\kappa = \langle \delta_k, \operatorname{Hess}[F] (\delta_k)_ {x} \rangle_{x}$,
+1. Set $\kappa = \langle \delta_k, \operatorname{Hess}[F] (\delta_k)_{x} \rangle_{x}$,
     $\alpha =\frac{\langle r_k, z_k \rangle_{x}}{\kappa}$ and
     $\langle \eta_k, \eta_k \rangle_{x}^{* } = \langle \eta_k, \operatorname{P}(\eta_k) \rangle_{x} +
     2\alpha \langle \eta_k, \operatorname{P}(\delta_k) \rangle_{x} +  {\alpha}^2
@@ -36,8 +36,8 @@ Repeat until a convergence criterion is reached
     return $\eta_{k+1} = \eta_k + \tau \delta_k$ and stop.
 3. Set $\eta_{k}^{* }= \eta_k + \alpha \delta_k$, if
     $\langle \eta_k, \eta_k \rangle_{x} + \frac{1}{2} \langle \eta_k,
-    \operatorname{Hess}[F] (\eta_k)_ {x} \rangle_{x} \leqq \langle \eta_k^{* },
-    \eta_k^{* } \rangle_{x} + \frac{1}{2} \langle \eta_k^{* },
+    \operatorname{Hess}[F] (\eta_k)_{x} \rangle_{x} \leqq \langle \eta_k^{* },
+    \eta_k^{*} \rangle_{x} + \frac{1}{2} \langle \eta_k^{* },
     \operatorname{Hess}[F] (\eta_k)_ {x} \rangle_{x}$
     set $\eta_{k+1} = \eta_k$ else set $\eta_{k+1} = \eta_{k}^{* }$.
 4. Set $r_{k+1} = r_k + \alpha \operatorname{Hess}[F] (\delta_k)_ {x}$,
@@ -51,6 +51,7 @@ Repeat until a convergence criterion is reached
 The result is given by the last computed $η_k$.
 
 ## Remarks
+
 The $\operatorname{P}(\cdot)$ denotes the symmetric, positive definite
 preconditioner. It is required if a randomized approach is used i.e. using
 a random tangent vector $\eta$ as initial
@@ -63,14 +64,16 @@ default, the preconditioner is just the identity.
 To step number 2: Obtain $\tau$ from the positive root of
 $\left\lVert \eta_k + \tau \delta_k \right\rVert_{\operatorname{P}, x} = \Delta$
 what becomes after the conversion of the equation to
+
 ````math
  \tau = \frac{-\langle \eta_k, \operatorname{P}(\delta_k) \rangle_{x} +
  \sqrt{\langle \eta_k, \operatorname{P}(\delta_k) \rangle_{x}^{2} +
  \langle \delta_k, \operatorname{P}(\delta_k) \rangle_{x} ( \Delta^2 -
  \langle \eta_k, \operatorname{P}(\eta_k) \rangle_{x})}}
  {\langle \delta_k, \operatorname{P}(\delta_k) \rangle_{x}}.
 ````
-It can occur that $\langle \delta_k, \operatorname{Hess}[F] (\delta_k)_ {x} \rangle_{x}
+
+It can occur that $\langle \delta_k, \operatorname{Hess}[F] (\delta_k)_{x} \rangle_{x}
 = \kappa \leqq 0$ at iteration $k$. In this case, the model is not strictly
 convex, and the stepsize $\alpha =\frac{\langle r_k, z_k \rangle_{x}}
 {\kappa}$ computed in step 1. does not give a reduction in the modelfunction
@@ -84,7 +87,7 @@ line. Thus when $\kappa \leqq 0$ at iteration k, we replace $\alpha =
 The other possibility is that $\eta_{k+1}$ would lie outside the trust-region at
 iteration k (i.e. $\langle \eta_k, \eta_k \rangle_{x}^{* }
 \geqq {\Delta}^2$ what can be identified with the norm of $\eta_{k+1}$). In
-particular, when $\operatorname{Hess}[F] (\cdot)_ {x}$ is positive definite
+particular, when $\operatorname{Hess}[F] (\cdot)_{x}$ is positive definite
 and $\eta_{k+1}$ lies outside the trust region, the solution to the
 trust-region problem must lie on the trust-region boundary. Thus, there
 is no reason to continue with the conjugate gradient iteration, as it
@@ -107,14 +110,8 @@ TruncatedConjugateGradientOptions
 ## Additional Stopping Criteria
 
 ```@docs
-stopIfResidualIsReducedByPower
-```
-```@docs
-stopIfResidualIsReducedByFactor
-```
-```@docs
+StopIfResidualIsReducedByPower
+StopIfResidualIsReducedByFactor
 StopWhenTrustRegionIsExceeded
-```
-```@docs
 StopWhenCurvatureIsNegative
 ```