The Riemannian Chambolle-Pock Algorithm
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
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
is fulfilled with)
- \[\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)\]
- \[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)\]
- 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}$
- \[\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.
Manopt.ChambollePock
— FunctionChambollePock(M, N, cost, x0, ξ0, m, n, prox_F, prox_G_dual, forward_operator, adjoint_DΛ)
Perform the Riemannian Chambolle–Pock algorithm.
Given a cost
function $\mathcal E\colon\mathcal M \to ℝ$ of the form
where $F\colon\mathcal M \to ℝ$, $G\colon\mathcal N \to ℝ$, and $\Lambda\colon\mathcal M \to \mathcal N$. The remaining input parameters are
x,ξ
primal and dual start points $x\in\mathcal M$ and $\xi\in T_n\mathcal N$m,n
base points on $\mathcal M$ and $\mathcal N$, respectively.forward_operator
the operator $Λ(⋅)$ or its linearization $DΛ(⋅)[⋅]$, depending on whether:exact
or:linearized
is chosen.adjoint_linearized_operator
the adjoint $DΛ^*$ of the linearized operator $DΛ(m)\colon T_{m}\mathcal M \to T_{Λ(m)}\mathcal N$prox_F, prox_G_Dual
the proximal maps of $F$ and $G^\ast_n$
By default, this performs the exact Riemannian Chambolle Pock algorithm, see the opional parameter DΛ
for ther linearized variant.
For more details on the algorithm, see[BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020].
Optional Parameters
acceleration
– (0.05
)dual_stepsize
– (1/sqrt(8)
) proximnal parameter of the primal proxΛ
(missing
) the exact operator, that is required if the forward operator is linearized;missing
indicates, that the forward operator is exact.primal_stepsize
– (1/sqrt(8)
) proximnal parameter of the dual proxrelaxation
– (1.
)relax
– (:primal
) whether to relax the primal or dualvariant
- (:exact
ifΛ
is missing, otherwise:linearized
) variant to use. Note that this changes the arguments theforward_operator
will be called.stopping_criterion
– (stopAtIteration(100)
) aStoppingCriterion
update_primal_base
– (missing
) function to updatem
(identity by default/missing)update_dual_base
– (missing
) function to updaten
(identity by default/missing)retraction_method
– (ExponentialRetraction()
) the rectraction to useinverse_retraction_method
- (LogarithmicInverseRetraction()
) an inverse retraction to use.vector_transport_method
- (ParallelTransport()
) a vector transport to use
Manopt.ChambollePock!
— FunctionChambollePock(M, N, cost, x, ξ, m, n, prox_F, prox_G_dual, forward_operator, adjoint_DΛ)
Perform the Riemannian Chambolle–Pock algorithm in place of x
, ξ
, and potenitally m
, n
if they are not fixed. See ChambollePock
for details and optional parameters.
Problem & Options
Manopt.PrimalDualProblem
— TypePrimalDualProblem {mT <: Manifold, nT <: Manifold} <: PrimalDualProblem} <: Problem
Describes a Problem for the linearized Chambolle-Pock algorithm.
Fields
M
,N
– two manifolds $\mathcal M$, $\mathcal N$cost
$F + G(Λ(⋅))$ to evaluate interims cost function valuesforward_oprator
the operator for the forward operation in the algorthm, either $Λ$ (exact) or $DΛ$ (linearized).linearized_adjoint_operator
The adjoint differential $(DΛ)^* \colon \mathcal N \to T\mathcal M$prox_F
the proximal map belonging to $f$prox_G_dual
the proximal map belonging to $g_n^*$Λ
– (fordward_operator
) for the linearized variant, this has to be set to the exact forward operator. This operator is required in several variants of the linearized algorithm. Since the exact variant is the default,Λ
is by default set toforward_operator
.
Constructor
LinearizedPrimalDualProblem(M, N, cost, prox_F, prox_G_dual, forward_operator, adjoint_linearized_operator,Λ=forward_operator)
Manopt.PrimalDualOptions
— TypePrimalDualOptions
A general type for all primal dual based options to be used within primal dual based algorithms
Manopt.ChambollePockOptions
— TypeChambollePockOptions <: PrimalDualOptions
stores all options and variables within a linearized or exact Chambolle Pock. The following list provides the order for the constructor, where the previous iterates are initialized automatically and values with a default may be left out.
m
- base point on $ \mathcal M $n
- base point on $ \mathcal N $x
- an initial point on $x^{(0)} \in \mathcal M$ (and its previous iterate)ξ
- an initial tangent vector $\xi^{(0)}\in T^*\mathcal N$ (and its previous iterate)xbar
- the relaxed iterate used in the next dual update step (when using:primal
relaxation)ξbar
- the relaxed iterate used in the next primal update step (when using:dual
relaxation)Θ
– factor to damp the helping $\tilde x$primal_stepsize
– (1/sqrt(8)
) proximal parameter of the primal proxdual_stepsize
– (1/sqrt(8)
) proximnal parameter of the dual proxacceleration
– (0.
) acceleration factor due to Chambolle & Pockrelaxation
– (1.
) relaxation in the primal relaxation step (to computexbar
)relax
– (_primal
) which variable to relax (:primal
or:dual
)stop
- aStoppingCriterion
type
– (exact
) whether to perform an:exact
or:linearized
Chambolle-Pockupdate_primal_base
((p,o,i) -> o.m
) function to update the primal baseupdate_dual_base
((p,o,i) -> o.n
) function to update the dual baseretraction_method
– (ExponentialRetraction()
) the rectraction to useinverse_retraction_method
- (LogarithmicInverseRetraction()
) an inverse retraction to use.vector_transport_method
- (ParallelTransport()
) a vector transport to use
where for the last two the functions a Problem
p
, Options
o
and the current iterate i
are the arguments. If you activate these to be different from the default identity, you have to provide p.Λ
for the algorithm to work (which might be missing
in the linearized case).
Constructor
ChambollePockOptions(m::P, n::Q, x::P, ξ::T, primal_stepsize::Float64, dual_stepsize::Float64;
acceleration::Float64 = 0.0,
relaxation::Float64 = 1.0,
relax::Symbol = :primal,
stopping_criterion::StoppingCriterion = StopAfterIteration(300),
variant::Symbol = :exact,
update_primal_base::Union{Function,Missing} = missing,
update_dual_base::Union{Function,Missing} = missing,
retraction_method = ExponentialRetraction(),
inverse_retraction_method = LogarithmicInverseRetraction(),
vector_transport_method = ParallelTransport(),
)
Useful Terms
Manopt.primal_residual
— Functionprimal_residual(p, o, x_old, ξ_old, n_old)
Compute the primal residual at current iterate $k$ given the necessary values $x_{k-1}, ξ_{k-1}, and $n_{k-1}$ from the previous iterate.
where $V_{\cdot\gets\cdot}$ is the vector transport used in the ChambollePockOptions
Manopt.dual_residual
— Functiondual_residual(p, o, x_old, ξ_old, n_old)
Compute the dual residual at current iterate $k$ given the necessary values $x_{k-1}, ξ_{k-1}, and $n_{k-1}$ from the previous iterate. The formula is slightly different depending on the o.variant
used:
For the :lineaized
it reads
and for the :exact
variant
where in both cases $V_{\cdot\gets\cdot}$ is the vector transport used in the ChambollePockOptions
.
Debug
Manopt.DebugDualBaseIterate
— FunctionDebugDualBaseIterate(io::IO=stdout)
Print the dual base variable by using DebugEntry
, see their constructors for detail. This method is further set display o.n
.
Manopt.DebugDualBaseChange
— FunctionDebugDualChange(a=StoreOptionsAction((:ξ)),io::IO=stdout)
Print the change of the dual base variable by using DebugEntryChange
, see their constructors for detail, on o.n
.
Manopt.DebugPrimalBaseIterate
— FunctionDebugPrimalBaseIterate(io::IO=stdout)
Print the primal base variable by using DebugEntry
, see their constructors for detail. This method is further set display o.m
.
Manopt.DebugPrimalBaseChange
— FunctionDebugPrimalBaseChange(a::StoreOptionsAction=StoreOptionsAction((:m)),io::IO=stdout)
Print the change of the primal base variable by using DebugEntryChange
, see their constructors for detail, on o.n
.
Manopt.DebugDualChange
— TypeDebugDualChange(opts...)
Print the change of the dual variable, similar to DebugChange
, see their constructors for detail, but with a different calculation of the change, since the dual variable lives in (possibly different) tangent spaces.
Manopt.DebugDualIterate
— FunctionDebugDualIterate(e)
Print the dual variable by using DebugEntry
, see their constructors for detail. This method is further set display o.ξ
.
Manopt.DebugDualResidual
— TypeDebugDualResidual <: DebugAction
A Debug action to print the dual residual. The constructor accepts a printing function and some (shared) storage, which should at least record :x
, :ξ
and :n
.
Manopt.DebugPrimalChange
— FunctionDebugPrimalChange(opts...)
Print the change of the primal variable by using DebugChange
, see their constructors for detail.
Manopt.DebugPrimalIterate
— FunctionDebugPrimalIterate(opts...)
Print the change of the primal variable by using DebugIterate
, see their constructors for detail.
Manopt.DebugPrimalResidual
— TypeDebugPrimalResidual <: DebugAction
A Debug action to print the primal residual. The constructor accepts a printing function and some (shared) storage, which should at least record :x
, :ξ
and :n
.
Manopt.DebugPrimalDualResidual
— TypeDebugPrimalDualResidual <: DebugAction
A Debug action to print the primaldual residual. The constructor accepts a printing function and some (shared) storage, which should at least record :x
, :ξ
and :n
.
Record
Manopt.RecordDualBaseIterate
— FunctionRecordDualBaseIterate(n)
Create an RecordAction
that records the dual base point, i.e. RecordEntry
of o.n
.
Manopt.RecordDualBaseChange
— FunctionRecordDualBaseChange(e)
Create an RecordAction
that records the dual base point change, i.e. RecordEntryChange
of o.n
with distance to the last value to store a value.
Manopt.RecordDualChange
— FunctionRecordDualChange()
Create the action either with a given (shared) Storage, which can be set to the values
Tuple, if that is provided).
Manopt.RecordDualIterate
— FunctionRecordDualIterate(ξ)
Create an RecordAction
that records the dual base point, i.e. RecordEntry
of o.ξ
, so .
Manopt.RecordPrimalBaseIterate
— FunctionRecordPrimalBaseIterate(x)
Create an RecordAction
that records the primal base point, i.e. RecordEntry
of o.m
.
Manopt.RecordPrimalBaseChange
— FunctionRecordPrimalBaseChange()
Create an RecordAction
that records the primal base point change, i.e. RecordEntryChange
of o.m
with distance to the last value to store a value.
Manopt.RecordPrimalChange
— FunctionRecordPrimalChange(a)
Create an RecordAction
that records the primal value change, i.e. RecordChange
, since we just redord the change of o.x
.
Manopt.RecordPrimalIterate
— FunctionRecordDualBaseIterate(x)
Create an RecordAction
that records the dual base point, i.e. RecordIterate
, i.e. o.x
.
Internals
Manopt.update_prox_parameters!
— Functionupdate_prox_parameters!(o)
update the prox parameters as described in Algorithm 2 of Chambolle, Pock, 2010, i.e.
- $\theta_{n} = \frac{1}{\sqrt{1+2\gamma\tau_n}}$
- $\tau_{n+1} = \theta_n\tau_n$
- $\sigma_{n+1} = \frac{\sigma_n}{\theta_n}$
- BergmannHerzogSilvaLouzeiroTenbrinckVidalNunez2020
R. Bergmann, R. Herzog, M. Silva Louzeiro, D. Tenbrinck, J. Vidal-Núñez: Fenchel Duality Theory and a Primal-Dual Algorithm on Riemannian Manifolds, arXiv: 1908.02022 accepted for publication in Foundations of Computational Mathematics
- 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