The Riemannian Chambolle-Pock algorithm

The Riemannian Chambolle-Pock algorithm is a generalization of the Chambolle-Pock algorithm Chambolle and Pock [CP11]. It is also known as primal-dual hybrid gradient (PDHG) or primal-dual proximal splitting (PDPS) algorithm.

In order to minimize over $p∈\mathcal M$ the cost function consisting of

\[F(p) + G(Λ(p)),\]

where $F:\mathcal M → \overline{ℝ}$, $G:\mathcal N → \overline{ℝ}$, and $Λ:\mathcal M →\mathcal N$. If the manifolds $\mathcal M$ or $\mathcal N$ are not Hadamard, it has to be considered locally only, that is on geodesically convex sets $\mathcal C \subset \mathcal M$ and $\mathcal D \subset\mathcal N$ such that $Λ(\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 Bergmann, Herzog, Silva Louzeiro, Tenbrinck and Vidal-Núñez [BHS+21]. In the following, the case of the exact, primal relaxed Riemannian Chambolle-Pock algorithm is described.

Given base points $m∈\mathcal C$, $n=Λ(m)∈\mathcal D$, initial primal and dual values $p^{(0)} ∈\mathcal C$, $ξ_n^{(0)} ∈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 algorithm performs the steps $k=0,1,…$ (until a StoppingCriterion is fulfilled)

\[ξ^{(k+1)}_n = \operatorname{prox}_{\tau_k G_n^*}\Bigl(ξ_n^{(k)} + \tau_k \bigl(\log_n Λ (\bar p^{(k)})\bigr)^\flat\Bigr)\]
\[p^{(k+1)} = \operatorname{prox}_{\sigma_k F}\biggl(\exp_{p^{(k)}}\Bigl( \mathcal P_{p^{(k)}\gets m}\bigl(-\sigma_k DΛ(m)^*[ξ_n^{(k+1)}]\bigr)^\sharp\Bigr)\biggr)\]
If $\gamma > 0$, update (otherwise keep $\theta_k = \theta_0$, $\sigma_{k+1} = \sigma_k$, $\tau_{k+1} = \tau_k$)
- $\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}$
\[\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 inverse 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 $Λ(m^{(k)})\neq n^{(k)}$ at some point. All these cases are covered in the algorithm.

Manopt.ChambollePock — Function

ChambollePock(M, N, f, p, X, m, n, prox_F, prox_G_dual, adjoint_linear_operator; kwargs...)
ChambollePock!(M, N, f, p, X, m, n, prox_F, prox_G_dual, adjoint_linear_operator; kwargs...)

Perform the Riemannian Chambolle-Pock algorithm.

Given a cost function $\mathcal{E}: \mathcal{M} → ℝ$ of the form

\[\mathcal{E}(p) = F(p) + G( Λ(p) ),\]

where $F:\mathcal{M} → ℝ$, $G:\mathcal{N} → ℝ$, and $Λ:\mathcal{M} → \mathcal{N}$.

This can be done inplace of $p$.

Input parameters

M::AbstractManifold: a Riemannian manifold $\mathcal{M}$
N::AbstractManifold: a Riemannian manifold $\mathcal{M}$
f: a cost function $f: \mathcal{M}→ ℝ$ implemented as (M, p) -> v
p::P: a point on the manifold $\mathcal{M}$
X::T: a tangent vector at the point $n$ on the manifold $\mathcal{N}$
m::P: a point on the manifold $\mathcal{M}$
n::P: a point on the manifold $\mathcal{N}$
adjoint_linearized_operator: the adjoint $DΛ^*$ of the linearized operator $DΛ: T_{m}\mathcal{M} → T_{Λ(m)}\mathcal{N}$
prox_F, prox_G_dual: the proximal maps of $F$ and $G^\ast_n$

If the forward operator Λ is provided, this performs the exact Riemannian Chambolle Pock algorithm; see the optional keyword linearized_forward_operator= for the linearized variant.

For more details on the algorithm, see [BHS+21].

Keyword Arguments

acceleration=0.05: acceleration parameter
callbacks::D = Dict{Symbol,Function}(): provided callback functions given either as a single function (symbol, problem, state, k) called in every hook or as a (vector of) pairs :hook => function, which are processed by process_callbacks_arg. As key you can either pass single symbol or an array of symbols to indicate a callback should be added in multiple places
dual_stepsize=1/sqrt(8): proximal parameter of the dual prox
evaluation::AbstractEvaluationType=AllocatingEvaluation(): specify whether the functions that return an array, for example a point or a tangent vector, work by allocating their result (AllocatingEvaluation) or whether they modify their input argument to return the result therein (InplaceEvaluation). Since usually the first argument is the manifold, the modified argument is the second.
inverse_retraction_method::AbstractInverseRetractionMethod=default_inverse_retraction_method(M, typeof(p)): an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inverses
inverse_retraction_method_dual::AbstractInverseRetractionMethod=default_inverse_retraction_method(N, typeof(n)): an inverse retraction $\operatorname{retr}^{-1}$ to use, see the section on retractions and their inverses
Λ=missing: the (forward) operator $Λ(⋅)$ (required for the :exact variant)
linearized_forward_operator=missing: its linearization $DΛ(⋅)[⋅]$ (required for the :linearized variant)
primal_stepsize=1/sqrt(8): proximal parameter of the primal prox
relaxation=1.: the relaxation parameter $θ$
relax=:primal: whether to relax the primal or dual
variant=:linearized if Λ is missing, otherwise :exact: variant to use. Note that this changes the arguments the forward_operator is called with.
stopping_criterion::StoppingCriterion=StopAfterIteration(200): a functor indicating that the stopping criterion is fulfilled
update_primal_base=missing: function to update m (identity by default/missing)
update_dual_base=missing: function to update n (identity by default/missing)
retraction_method::AbstractRetractionMethod=default_retraction_method(M, typeof(p)): a retraction $\operatorname{retr}$ to use, see the section on retractions
vector_transport_method::AbstractVectorTransportMethod=default_vector_transport_method(M, typeof(p)): a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports
vector_transport_method_dual::AbstractVectorTransportMethod=default_vector_transport_method(N, typeof(n)): a vector transport $\mathcal T_{⋅←⋅}$ to use, see the section on vector transports

Output

The obtained approximate minimizer $p^*$. To obtain the whole final state of the solver, see get_solver_return for details, especially the return_state= keyword.